In[]:=
CompoundExpression[
]
​​deploy
Mon 18 Dec 2023 13:56:28
Parent pseudospectra.nb
​
Lyapunov Theory for Discrete Time Systems notability
​
​

Util

In[]:=
exportAnim[fn_,images_,is_]:=Module[{},​​rasters=Rasterize[Style[#,AntialiasingTrue],ImageSizeis,RasterSize2is]&/@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
d_
A_
⟺
distinguishMatrixPowerAndPower[A_,d_]
;​​*)​​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[
-1
X
];​​DA=DiagonalMatrix[Eigenvalues@X];​​DB=DiagonalMatrix[Eigenvalues[
-1
X
]];​​-T.
-1
T
.(
-1
X
).U
DA.ones-ones.DB
.
-1
U
//Chop​​;​​​​(* factor=1: identity transformation, factor=infinity -- approaching divergent *)​​makeRapetti[factor_]:=​​ii=IdentityMatrix[2];​​PA={{2, 100},{-300,100}};​​alpha=Min
2Re[#]
2
Abs[#]
&/@Eigenvalues[PA];​​Nii-1-
1
factor
alpha PA​​;​​​​​​lyapunovCondition[mat_]:=Module[{P,ii},​​ii=IdentityMatrix[First@Dimensions@mat];​​P=DiscreteLyapunovSolve[mat,-ii];​​LinearAlgebra`Private`MatrixConditionNumber[P]​​];

When is (I-x A) contractive?


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 eigenvalues change

Take Rapetti example, make it increasingly unstable, track eigenvalues,singular values

How norms/spectra change with preconditioning

Old Solving for coordinate change/grid

Off-diagonal example

How things diverge