Solving Autonomous Linear Systems

Eigenvalues and eigenvectors

The simplest way to solve a linear system of ODE’s is to find n linearly independent eigenvectors.  This is generally possible, in principle, but finding eigenvalues and eigenvectors is complicated.

Here is a simple example.  We enter a matrix and find eigenvalues and eigenvectors.

A= {{7,-8,-2},{10,-9,-4},{0,-4,1}};A//MatrixForm
CharacteristicPolynomial[A,z]
%//Factor
evals=Eigenvalues[A]
evects=Eigenvectors[A]

Two of the eigenvalues are complex as are the corresponding eigenvectors.  The general complex solution would be an arbitrary complex linear combination of these.  But it is hard to see the form of the general real solution from this.

complexsol = c1*Exp[t*evals[[1]]]*evects[[1]]+ c2*Exp[t*evals[[2]]]*evects[[2]]+ c3*Exp[t*evals[[3]]]*evects[[3]]

We can get the general real form by constructing a principal fundamental matrix.  First make a complex fundamental matrix out of the eigenvalues and eigenvectors, then make a principal one from that.

P = evects//Transpose;P//MatrixForm
phi = P.DiagonalMatrix[Exp[t*evals]];phi//MatrixForm
psi = P.DiagonalMatrix[Exp[t*evals]].Inverse[P]//ComplexExpand//Simplify

The principal fundamental matrix of an autonomous system is actually the matrix exponential.  Mathematica even has a command for this.

MatrixExp[t*A]

Diagonalization and Real Canonical Form

Another point of view about this approach is to use the eigenvalues and eigenvectors to transform the differential equation into a simpler form.  Using complex vectors we can usually diagonalize the matrix.

B= Inverse[P].A.P//Simplify;
B//MatrixForm

Alternatively if we use the real and imaginary parts of the complex pairs of eigenvectors we can get a simple real matrix.  From this real canonical form we can easily understand the flow in the new coordinates.  The flow in the old coordinates is related by the linear transformation with matrix Q.

u=Re[evects[[2]]];
v=Im[evects[[2]]];
Q= {evects[[1]],u,v}//Transpose;
Q//MatrixForm
B1=Inverse[Q].A.Q//Simplify;
B1//MatrixForm

Phase Portraits

Here is one way to make a phase portait.  Just compute several well-chosen orbits and plot using ParametricPlot or ParametricPlot3D.    First plot the canonical form ODE.

ode= MapThread[Equal,{{x1'[t],x2'[t],x3'[t]},B1.{x1[t],x2[t],x3[t]}}]
ic[{x10_,x20_,x30_}] := {x1[0]==x10,x2[0]==x20,x3[0]==x30}

Choose some initial conditions in the invariant eigensubpaces and others not in these spaces.

sol1= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[0.01*{0,1,0}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol2= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[0.01*{0,0,1}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol3= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[-0.01*{0,1,0}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol4= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[-0.01*{0,0,1}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol5= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[5*{1,0,0}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol6= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[-5*{1,0,0}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol7= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[5*{1,0,0}+{0,1,0}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];
sol8= {x1[t],x2[t],x3[t]}/.NDSolve[ode~Join~ic[-5*{1,0,0}+{0,0,1}],{x1[t],x2[t],x3[t]},{t,0,10}][[1]];

ParametricPlot3D[{sol1,sol2,sol3,sol4,sol5,sol6,sol7,sol8},{t,0,10},PlotRange->{{-5,5},{-5,5},{-5,5}}]

The solutions of the original system are the images of these under the linear map Q.

ParametricPlot3D[{Q.sol1,Q.sol2,Q.sol3,Q.sol4,Q.sol5,Q.sol6,Q.sol7,Q.sol8},{t,0,10},PlotRange->{{-5,5},{-5,5},{-5,5}}]