In[]:=
Fri 24 Nov 2023 23:57:44
Rapetti example
Rapetti example
In[]:=
ii=IdentityMatrix[2];PA={{2,100},{-300,100}};rho[mat_]:=Max[Abs/@Eigenvalues[mat]];alpha=Min&/@Eigenvalues[PA]rho[ii-alphaPA]
2Re[#]
2
Abs[#]
Out[]=
51
15100
Out[]=
1
In[]:=
T=ii-alphaPA;numSteps=200;numPoints=20;stepLength=1/20;xs0=CirclePoints[numPoints]//N;xs=Table[Transpose[Re@MatrixPower[T,k].Transpose[xs0]],{k,0,numSteps*stepLength,stepLength}];Print[xs//Dimensions]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]];combinedTrajectoryPlot[s_]:=Show[contour,ListPlot[xs[[s]]],drawTrajectories[s],PlotRange->{{-1.3,1.3},{-1.3,1.3}},AspectRatio->1,PlotLabel->StringForm["values of x",MatrixForm[nf[T]],NumberForm[N[s*stepLength],{3,2}]]];plots=Table[combinedTrajectoryPlot[s],{s,1,Length[xs]}];ListAnimate[plots]
``
``
{201,20,2}
Out[]=
Remove imaginary component
Remove imaginary component
In[]:=
T=ii-alphaPA;{evals,evecs}=Eigensystem[T];U=Transpose[evecs];(*Norm[T-U.DiagonalMatrix[evals].Inverse[U]]*)T2=U.DiagonalMatrix[Re@evals].Inverse[U];T=T2;numSteps=200;numPoints=20;stepLength=1/2;xs0=CirclePoints[numPoints]//N;xs=Table[Transpose[Re@MatrixPower[T,k].Transpose[xs0]],{k,0,numSteps*stepLength,stepLength}];Print[xs//Dimensions]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]];combinedTrajectoryPlot[s_]:=Show[contour,ListPlot[xs[[s]]],drawTrajectories[s],PlotRange->{{-1.3,1.3},{-1.3,1.3}},AspectRatio->1,PlotLabel->StringForm["values of x",MatrixForm[nf[T]],NumberForm[N[s*stepLength],{3,2}]]];plots=Table[combinedTrajectoryPlot[s],{s,1,Length[xs]}];ListAnimate[plots]
``
``
{201,20,2}
Out[]=
Remove real component
Remove real component
T=ii-alphaPA;{evals,evecs}=Eigensystem[T];U=Transpose[evecs];(*Norm[T-U.DiagonalMatrix[evals].Inverse[U]]*)T2=U.DiagonalMatrix[evals-Re[evals]].Inverse[U];T=N@T2;Print["T2=",T2//N//Chop//MatrixForm];nf[x_]:=NumberForm[Chop@N@x,{3,3}];numSteps=200;numPoints=20;stepLength=1/2;xs0=CirclePoints[numPoints]//N;xs=Table[Transpose[Re@MatrixPower[T,k].Transpose[xs0]],{k,0,numSteps*stepLength,stepLength}];Print[xs//Dimensions]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]];combinedTrajectoryPlot[s_]:=Show[contour,ListPlot[xs[[s]]],drawTrajectories[s],PlotRange->{{-1.3,1.3},{-1.3,1.3}},AspectRatio->1,PlotLabel->StringForm["values of x",MatrixForm[nf[T]],NumberForm[N[s*stepLength],{3,2}]]];plots=Table[combinedTrajectoryPlot[s],{s,1,Length[xs]}];ListAnimate[plots]
``
``
T2=
0.165497 | -0.337748 |
1.01325 | -0.165497 |
{201,20,2}
Out[]=
Amplify imaginary component
Amplify imaginary component