In[]:=
deploy
Tue 12 Jul 2022 17:02:49
Discrete Lyapunov equation
Discrete Lyapunov equation
In[]:=
On[Assert];vec[W_]:=Transpose@{Flatten@Transpose[W]};unvec[Wf_, rows_]:=Transpose[Flatten/@Partition[Wf,rows]];unvec[Wf_]:=Module[{d},d=Sqrt[Dimensions[Wf]//First];Assert[dFloor[d]];Assert[d*dFirst@Dimensions[Wf]];unvec[Wf,d]];CircleTimes=KroneckerProduct;Needs["Notation`"];distinguishMatrixPowerAndPower[A_,d_]:=If[SquareMatrixQ[A],MatrixPower[A,d],Power[A,d]];Notation ⟺ ;rvec[mat_]:=Transpose[Eigenvectors@mat]; (* right eigenvectors as columns *)lvec[mat_]:=Transpose[Eigenvectors[mat]]; (* left eigenvectors as columns *)
Least squares form
Least squares form
In[]:=
solveStein[A_,B_,E_]:=Module[{},ii=IdentityMatrix[First[Dimensions[A]]];X=unvec@LeastSquares[B⊗A-ii⊗ii,vec[E]];Print["errorLsqr=",Norm[A.X.B-X-E]];X];n=2;A0={{-1,4},{3,0}};B0={{-4,3},{2,6}};C0={{-5,0},{-2,-6}};DiscreteLyapunovSolve[A0,B0,C0]solveStein[A0,B0,C0]
Out[]=
-,-,,-
903
43175
14469
43175
10372
43175
9519
43175
errorLsqr=0
Out[]=
-,-,,-
903
43175
14469
43175
10372
43175
9519
43175
Eigenvalue form
Eigenvalue form
In[]:=
A0={{-1,4},{3,0}};B0={{-4,3},{2,6}};C0={{-5,0},{-2,-6}};ones=ConstantArray[1,{Length[A0],Length[A0]}];Avals=Eigenvalues[A0];Bvals=Eigenvalues[];Avec=rvec[A0];Bvec=rvec[];DA=DiagonalMatrix[Avals];DB=DiagonalMatrix[Bvals];sol0=Avec..C0.Inverse[B0].Bvec.;Print[sol0//Simplify];sol1=DiscreteLyapunovSolve[A0,B0,C0];Print[sol1//Simplify];Reduce[sol0==sol1]
-1
B0
-1
B0
-1
Avec
DA.ones-ones.DB
-1
Bvec
-,-,,-
903
43175
14469
43175
10372
43175
9519
43175
-,-,,-
903
43175
14469
43175
10372
43175
9519
43175
Out[]=
True