In[]:=
Mon 18 Dec 2023 13:56:28
Parent pseudospectra.nb
Lyapunov Theory for Discrete Time Systems notability
Lyapunov Theory for Discrete Time Systems notability
Util
Util
In[]:=
exportAnim[fn_,images_,is_]:=Module[{},rasters=Rasterize[Style[#,AntialiasingTrue],ImageSizeis,RasterSize2is]&/@images;SetDirectory[NotebookDirectory[]];SetDirectory["export"];Export[fn,rasters,"AnimationRepetitions"Infinity]];(*Needs["Notation`"];distinguishMatrixPowerAndPower[A_,d_]:=If[SquareMatrixQ[A],MatrixPower[A,d],Power[A,d]];Notation ⟺ ;*)LSolve[X_]:=Module{rvec,lvec,normlyap,n,ones,T,U,DA,DB},rvec[mat_]:=Transpose[Eigenvectors@mat]; (* right eigenvectors as columns *)lvec[mat_]:=Transpose[Eigenvectors[mat]]; (* left eigenvectors as columns *)(* approximate equality testing *)DotEqual[a_,b_]:=Norm[Flatten[{a}]-Flatten[{b}],∞]<1*^-9;normlyap[mat_]:=With[{x=mat/Norm[mat,"Frobenius"]},x/Sign[First@First@x]];n=First@Dimensions[X];ones=Array[1&,{n,n}];T=rvec[X];U=lvec[];DA=DiagonalMatrix[Eigenvalues@X];DB=DiagonalMatrix[Eigenvalues[]];-T..().U.//Chop;(* factor=1: identity transformation, factor=infinity -- approaching divergent *)makeRapetti[factor_]:=ii=IdentityMatrix[2];PA={{2, 100},{-300,100}};alpha=Min&/@Eigenvalues[PA];Nii-1-alpha PA;lyapunovCondition[mat_]:=Module[{P,ii},ii=IdentityMatrix[First@Dimensions@mat];P=DiscreteLyapunovSolve[mat,-ii];LinearAlgebra`Private`MatrixConditionNumber[P]];
-1
X
-1
X
-1
T
-1
X
DA.ones-ones.DB
-1
U
2Re[#]
2
Abs[#]
1
factor
When is (I-x A) contractive?
When is (I-x A) contractive?
Lyapunov contractive solve of Rapeti with grid
Lyapunov contractive solve of Rapeti with grid
Summary: Lyapunov preconditioning transforms it into convergent matrix, but it’s still not normal, so rates of convergence for max/min forever oscillate between eigenvalues. Rates of convergence for median/mean are close to spectral radius.
ClearAll["Global`*"];
In[]:=
(*Symmetricsquareroot*)symsqrt[m_]:=Module[{U,S,W},{U,S,W}=SingularValueDecomposition[m];U.Sqrt[S].Transpose[W]];nf[num_]:=NumberForm[N@num,{3,2}];(*setfactorto∞togetnon-convergent*)ii=IdentityMatrix[2]//N;T=makeRapetti[2];(*Numberofstepstorunanimationfor*)numSteps=10;numPoints=5;(*Howmuchtoincrementmatrixpowerateachstep*)stepLength=1/2;xs0=CirclePoints[numPoints]//N;xs=Table[Transpose[Re@MatrixPower[T,k].Transpose[xs0]],{k,0,numSteps*stepLength,stepLength}];xs=Table[xs0.MatrixPower[T,k],{k,0.,numSteps,stepLength}];(*stepsxbatchxdims*)(*xs=Table[xs0.MatrixPower[T,k],{k,1.,20,1/2}];*)drawTrajectory[s_]:=Graphics[{Opacity[.1],Point[xs[[;;s]]]}];drawTrajectories[s_]:=Graphics[{Opacity[.1]}~Join~(Line/@Transpose[xs][[All,;;s,All]])];contour=ContourPlot[x^2+y^2==1,{x,-1,1},{y,-1,1},ContourStyle->Directive[Opacity[.4],Gray]];P=DiscreteLyapunovSolve[T,-ii];sqrtP=CholeskyDecomposition[P];(*whatifweuseadiffone?*)sqrtP=symsqrt[P];isqrtP=Inverse[symsqrt@P];(*Plotbound*)bound=1;(*creategridforforgivenparams*)gridBound=2;gridWidth=0.1;gridTransform=sqrtP;xgrid:=ParametricPlot@@{Table[gridTransform.{x,yval},{yval,Range[-gridBound,gridBound,gridWidth]}],{x,-3,3},PlotStyle->Directive[Blue,Opacity[.3]],PlotRange->{{-bound,bound},{-bound,bound}}};ygrid:=ParametricPlot@@{Table[gridTransform.{xval,y},{xval,Range[-gridBound,gridBound,gridWidth]}],{y,-3,3},PlotStyle->Directive[Blue,Opacity[.3]],PlotRange->{{-bound,bound},{-bound,bound}}};fullGrid:=Show[xgrid,ygrid,PlotRange->{{-bound,bound},{-bound,bound}}];combinedTrajectoryPlot[s_]:=Show[fullGrid,contour,ListPlot[xs[[s]],PlotStyle->Red],drawTrajectories[s],PlotRange->{{-bound,bound},{-bound,bound}},AspectRatio->1,Axes->None,PlotLabel->StringForm["values of x",MatrixForm[nf@T],nf[s*stepLength]]];ListLinePlot[Map[#.Inverse[P].#&,xs,{2}],PlotLabel->"Energies",AxesLabel->{"time","energy"}]plots1=Table[combinedTrajectoryPlot[s],{s,1,Length[xs]}];ListAnimate[plots1]bound=Norm[isqrtP];T=Inverse[sqrtP].T.sqrtP;(*ys=Map[Pisqrt.#&,xs,{2}];*)xs0=xs0.Inverse[sqrtP]//N;xs=Table[Transpose[Re@MatrixPower[T,k].Transpose[xs0]],{k,0,numSteps*stepLength,stepLength}];xs=Table[xs0.MatrixPower[T,k],{k,0.,numSteps,stepLength}];contour=ContourPlot[{x,y}.P.{x,y}==1,{x,-bound,bound},{y,-bound,bound},ContourStyle->Directive[Opacity[.4],Gray]];gridBound=2;gridWidth=0.1;gridTransform=IdentityMatrix[2];fullGrid=Show[xgrid,ygrid,PlotRange->{{-bound,bound},{-bound,bound}}];plots2=Table[combinedTrajectoryPlot[s],{s,1,Length[xs]}];ListAnimate[plots2]
``
``
How Lyapunov preconditioned matrix converges
How Lyapunov preconditioned matrix converges
How eigenvalues change
How eigenvalues change
Take Rapetti example, make it increasingly unstable, track eigenvalues,singular values
How norms/spectra change with preconditioning
How norms/spectra change with preconditioning
Old Solving for coordinate change/grid
Old Solving for coordinate change/grid
Off-diagonal example
Off-diagonal example
How things diverge
How things diverge