Linear Systems

Transforming a Sysytem into an N-th order Equation

Solve the initial value problem
    x1'[t] ==  5 x1[t] + 3 x2[t],
    x2'[t] == -3 x1[t] - 2 x2[t],
    x1[0] == 0, x2[0] == 1.
by transforming it into a second order equation as in a) and plot the solution for t>=0.

Solution

Let y[t] == x1[t], y'[t] == x1'[t]. From the first equation:

x2[t]  == (x1'[t] - 5x1[t])/3 == (y'[t] - 5y[t])/3;
x2'[t] == (y''[t] - 5y'[t])/3

From the second equation,

x2'[t] == -3x1[t] - 2x2[t] == -3y[t] - 2(y'[t] - 5y[t])/3

Setting the two equal gives

(y''[t] - 5y'[t])/3 == -3y[t] - 2(y'[t] - 5y[t])/3
y''[t] - 3 y'[t] - y[t] == 0

The initial conditions for y[t] are: y[0] == x1[0] == 0; from the first equation:
y'[0] == 5x1[0] + 3x2[0] == 3.

Clear[y,t]

soln = DSolve[{y''[t] - 3 y'[t] - y[t] == 0,
               y[0] == 0, y'[0] == 3}, y[t],t]

y[t_] = y[t] /.First[soln]

The solution to the original system is:

x1[t_] = y[t]
x2[t_] = (y'[t] - 5y[t])/3 // Simplify

ParametricPlot[{x1[t],x2[t]}, {t,0,1}]

In one step:

Clear[x1,x2,t]

DSolve[{x1'[t] ==  5 x1[t] + 3 x2[t],
        x2'[t] == -3 x1[t] - 2 x2[t],
        x1[0] == 0, x2[0] == 1}, {x1[t],x2[t]}, t]//Simplify

Solving Systems of Equations

Solve the system of linear equations
    3x1 + 2x2 +  x3       ==  4,
    2x1 +  x2       +  x4 ==  4,
     x1 +        x3 + 2x4 ==  4,
           x2 + 2x3 + 3x4 ==  4,
by using Solve[], LinearSolve[], and by using the inverse of the matrix of coefficients.

Solution

Solve[{3x1 + 2x2 +  x3       ==  4,
       2x1 +  x2       +  x4 ==  4,
        x1 +        x3 + 2x4 ==  4,
              x2 + 2x3 + 3x4 ==  4}]

?LinearSolve

LinearSolve[{{3,2,1,0},
            {2,1,0,1},
            {1,0,1,2},
            {0,1,2,3}},{4,3,2,1}]

The matrix of coefficients is

A = {{3,2,1,0},
     {2,1,0,1},
     {1,0,1,2},
     {0,1,2,3}};

The inverse of A is:

B = Inverse[A];
B//MatrixForm

Applying this inverse to the vector on the RHS of the equation gives the same solution as LinearSolve[]:

B.{4,3,2,1}

Finding Eigenvalues and Eigenvectors

Find the eigenvalues and eigenvectors of the matrix
    A = {{4, 3, -1},
         {2, 5, -1},
         {4, 6, 0}}
Compare your results to the output of Eigensystem[]. Then find a matrix T such that Inverse[T].A.T == d where d is a diagonal matrix.

Solution

A = {{4, 3, -1},
      {2, 5, -1},
      {4, 6, 0}};

First compute the eigenvalues by finding the roots of

Det[A - r*IdentityMatrix[3]] == 0

Solve[Det[A - r*IdentityMatrix[3]] == 0]

Next find the eigenvectors corresponding to each of the eigenvalues by solving the equations.

eqns = (A-r*IdentityMatrix[3]).{v1,v2,v3} == {0,0,0} //Simplify

For the eigenvalue 2:

Solve[eqns /. {r->2}]

The values of v2 and v3 can be chosen freely, say constants 2a and 2b, respectively. The solution vector is then
    {-3a + b, 2a, 2b} == a{-3,2,0} + b{1,0,2}
Two eigenvectors corresponding to the repeated eigenvalue 2 are thus

X1 = {-3,2,0};  X2 = {1,0,2};

For the eigenvalue 5:

Solve[eqns /. {r->5}]

The value of v3 can be chosen freely, say 2c. The solution vector is then
    {c,c,2c} == c{1,1,2}
An eigenvector corresponding to the eigenvalue 5 is thus

X3 = {1,1,2};

Let's check that the vectors we have found are in fact eigenvectors:

A.X1 == 2 X1

A.X2 == 2 X2

A.X3 == 5 X3

Finding eigenvalues and eigenvectors can be done in one step:

?Eigensystem

Eigensystem[A]

By chance Eigensystem[] gives the same two eigenvectors for the repeated eigenvalue 2. There are infintely many choices possible for these two vectors by varying the constants a and b above. For example:

X11 = 2{-3,2,0} + {1,0,2}

X22 = {-3,2,0} + 2{1,0,2}

A.X11 == 2 X11

A.X22 == 2 X22

To get a diagonal matrix d == Inverse[T].A.T, just let T be the matrix whose columns are the eigenvectors X1, X2, and X3, see p. 362.

T = Transpose[{X1,X2,X3}];
T//MatrixForm

d = Inverse[T].A.T;
d//MatrixForm

The diagonal entries are the corresponding eigenvectors.

Wronskians and Solutions of Systems of Differential Equations

Calculate the Wronskian of the vector functions
    X1 = {t^2-t, t, t-1};
    X2 = {  0,  t^3, t };
    X3 = {  t,  t^2, 0 };
For what values of t are the vectors x1, x2, x3 linearly independent? On what intervals are the vector functions x1, x2, x3 linearly independent? Find a homogeneous system of first order equations satisfied by x1, x2, x3 and check where the coefficients are continuous.

Solution

First we caclulate the Wronskian of the vector functions.

Clear[X1,X2,X3]

X1 = {t^2-t, t, t-1};
X2 = {  0,  t^3, t };
X3 = {  t,  t^2, 0 };

w = Det[{X1,X2,X3}]

Solve[w == 0]

As vectors, X1, X2, X3 are linearly independent as long as t is not 0, or
(2 ± 2*3^(1/2))/4. As vector functions they are linearly independent  on any interval. The only linear comination that equals {0,0,0} for all t is when the coeficients are all zero:

eqns = { c1 X1 + c2 X2 + c3 X3 == {0,0,0} /. t->0,
         c1 X1 + c2 X2 + c3 X3 == {0,0,0} /. t->1,
         c1 X1 + c2 X2 + c3 X3 == {0,0,0} /. t->2}

Solve[eqns]

To find a homogeneous system that is satisfied by X1, X2, X3, we need only solve the following equations for the coefficients:

P = {{p11,p12,p13},
     {p21,p22,p23},
     {p31,p32,p33}};

eqns = { D[X1,t] == P.X1,
         D[X2,t] == P.X2,
         D[X3,t] == P.X3}

Solve[eqns]//Simplify

We see that the coefficients are not all continuous at t == 0 and (1 ± 3^(1/2))/2. Note that if they were continuous at these points, then Theorem 7.4.3 would imply that the Wronskian would not vanish for these values of t.

Solving Homogeneous Systems with Constant Coefficients and Plotting Trajectories

Saddle (eigenvalues of different sign)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{ 7 , 5},
     {-10,-8}};

Solution

Find the eigenvalues and eigenvectors of A:

Eigensystem[A]

This gives the two fundamental solutions

X1[t_] = {-1,1}E^(2t);
X2[t_] = {-1,2}E^(-3t);

The general solution is thus

X[t_] = c1 X1[t] + c2 X2[t]

Let us plot the fundamental solutions, whose trajectories are lines, in different colors. Using Evaluate[] keeps Mathematica from complaining about not being able to compile the function.

p1 = ParametricPlot[Evaluate[{X1[t],-X1[t]}],{t,-2,2},
     PlotStyle->RGBColor[1,0,0]]

p2 = ParametricPlot[Evaluate[{X2[t],-X2[t]}],{t,-2,2},
     PlotStyle->RGBColor[0,0,1]]

We can use Table[] to generate a list of solutions for various choices of c1 and c2. We must use Flatten[...,1]to remove an  extra layer of braces on the output of Table[].

solns = Flatten[Table[X[t],{c1,-2,2}, {c2,-2,2}],1];

p = ParametricPlot[Evaluate[solns], {t,-2,2},
    PlotRange->{{-5,5},{-5,5}}]

Now superimpose the fundamental solutions.

Show[p,p1,p2, PlotRange->{{-5,5},{-5,5}}]

The plots show that almost all solutions approach infinity asymtotic to the solution X1[t] (the red line). Solutions that start on the blue line (multiples of X2[t]) are the only exceptions and they approach the origin.

Note that the system can also be solved using DSolve[], although the constants come out arranged differently than above.

Clear[x1,x2,x3]

DSolve[{x1'[t] ==  7x1[t] + 5x2[t],
        x2'[t] ==-10x1[t] - 8x2[t]}, {x1[t],x2[t]}, t]

Stable Node (eigenvalues both negative)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{1,-2},
     {3,-4}};

Solution

We follow the same steps as in the first example:

Eigensystem[A]

This gives the two fundamental solutions

X1[t_] = {1,1}E^(-t);
X2[t_] = {2,3}E^(-2t);

The general solution is thus

X[t_] = c1 X1[t] + c2 X2[t]

Let us plot the fundamental solutions in red and blue as before.

p1 = ParametricPlot[Evaluate[{X1[t],-X1[t]}],{t,-1,3},
     PlotStyle->RGBColor[1,0,0]]

p2 = ParametricPlot[Evaluate[{X2[t],-X2[t]}],{t,-1,3},
     PlotStyle->RGBColor[0,0,1]]

solns = Flatten[Table[X[t],{c1,-2,2}, {c2,-2,2}],1];

We'll keep the plot within -1 and 1 in both directions so we can better see what is happening at the origin.

p = ParametricPlot[Evaluate[solns], {t,-1,3},
PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic]

Now superimpose the fundamental solutions.

Show[p,p1,p2, PlotRange->{{-1,1},{-1,1}}]

From these graphs we see that all solutions approach the origin tangent to the fundamental solution X1[t] (the red line).

Unstable Node (eigenvalues both positive)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{5,-1},
     {3, 1}};

Then find and plot the solution satisfyting the initial condition X[0] = {2,-1}.

Solution

Eigensystem[A]

This gives the two fundamental solutions

X1[t_] = {1,1}E^(4t);
X2[t_] = {1,3}E^(2t);

The general solution is

X[t_] = c1 X1[t] + c2 X2[t]

p1 = ParametricPlot[Evaluate[{X1[t],-X1[t]}],{t,-3,3},
     PlotStyle->RGBColor[1,0,0]]

p2 = ParametricPlot[Evaluate[{X2[t],-X2[t]}],{t,-6,6},
     PlotStyle->RGBColor[0,0,1]]

solns = Flatten[Table[X[t],{c1,-3,3},{c2,-3,3}],1];

p = ParametricPlot[Evaluate[solns],{t,-3,3},AspectRatio->1]

Now superimpose the fundamental solutions. Since the soutions go to infinity, we use a large plot range to see their behaviour for large t.

Show[p,p1,p2, PlotRange->{{-100,100},{-100,100}},
AspectRatio->1]

From these graphs we see that all solutions approach infinity parallel to the solution X1[t] (the red line). The solution satisfying the given initial condition is found by solving:

X[0] == {2,-1}

Solve[X[0] == {2,-1}]

X3[t_] = X[t] /. First[%]

p3=ParametricPlot[X3[t],{t,0,3},PlotStyle->RGBColor[0,1,0]]

Show[p1,p2,p3,PlotRange->{{-100,100},{-100,100}},AspectRatio->1]

Non-isolated fixed points (zero eigenvalue)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{ 2,10},
     {-1,-5}};

Solution

Eigensystem[A]

X1[t_] = {-5,1}E^(0*t)

X2[t_] = {-2,1}E^(-3t);

The general solution is

X[t_] = c1 X1[t] + c2 X2[t]

The trajectory of X1[t] is a fixed point. The different multiples of X1[t] give a line of fixed points. The trajectory X2[t] move along the line y = -(1/2)x towards the origin.

p1 = ParametricPlot[Evaluate[{c1 X1[t]}],{c1,-5,5},
     PlotStyle->RGBColor[1,0,0]]

p2 = ParametricPlot[Evaluate[{X2[t],-X2[t]}],{t,-3,3},
     PlotStyle->RGBColor[0,0,1]]

solns = Flatten[Table[X[t],{c1,-4,4},{c2,-1,1}],1];

p = ParametricPlot[Evaluate[solns], {t,-3,3}]

Show[p,p1,p2,PlotRange->{{-20,20},{-10,10}}]

All the non-fixed point trajectories move along lines parallel to the trajectory X2[t] (the blue line) and approach a fixed point (i.e. a multiple of X2[t] on the red line).

Center (imaginary eigenvalues)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{ 1, 5},
     {-2,-1}};

Solution

Eigensystem[A]

This gives two conjugate complex-valued solutions. The real and imaginary part of either one then gives two linearly independent real-valued solutions.

Z[t_] = {-(1/2)-(3/2)I, 1}E^(3I*t);

The real and imaginary parts are

X1[t_] = Re[Z[t]]//ComplexExpand

X2[t_] = Im[Z[t]]//ComplexExpand

The general solution is

X[t_] = c1 X1[t] + c2 X2[t]

The trajectories of X1[t] and X2[t] are the same, both ellipses, even though they are linearly independent (they trace the path at different times).

p1 = ParametricPlot[Evaluate[X1[t]],{t,-Pi/3,Pi/3},
     AspectRatio->Automatic]

p2 = ParametricPlot[Evaluate[X2[t]],{t,-Pi/3,Pi/3},
     AspectRatio->Automatic]

solns = Table[X[t]/.c2->0,{c1,1,4}];

p = ParametricPlot[Evaluate[solns],{t,-Pi/3,Pi/3},
    AspectRatio->Automatic]

Stable Spiral (complex eigenvalues, negative real part)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{-1,-4},
     { 1,-1}};

Solution

Eigensystem[A]

This gives two conjugate complex-valued solutions. The real and imaginary part of either one then gives two linearly independent real-valued solutions.

Z[t_] = {2I, 1}E^((-1+2I)*t);

The real and imaginary parts are

X1[t_] = Re[Z[t]]//ComplexExpand

X2[t_] = Im[Z[t]]//ComplexExpand

The general solution is

X[t_] = c1 X1[t] + c2 X2[t]

The trajectories of X1[t] and X2[t] spiral into the origin.

p1 = ParametricPlot[Evaluate[X1[t]],{t,-3,3}]

p1 = ParametricPlot[Evaluate[X2[t]],{t,-3,3}]

solns = Flatten[Table[X[t],{c1,-2,2},{c2,-2,2}],1];

p = ParametricPlot[Evaluate[solns],{t,-3,3},
    PlotRange->{{-3,3},{-3,3}},AspectRatio->Automatic]

The plots show that all solutions spiral into the origin.

Unstable Spiral (complex eigenvalues, positive real part)

Find the general solution and plot several trajectories of X'[t] = A.X[t] where

A = {{ 3,-4},
     { 2,-1}};

Solution

Eigensystem[A]

A complex-valued solution is

Z[t_] = {1 + I, 1}E^((1 + 2I)*t);

The real and imaginary parts are

X1[t_] = Re[Z[t]]//ComplexExpand

X2[t_] = Im[Z[t]]//ComplexExpand

The general solution is

X[t_] = c1 X1[t] + c2 X2[t]

The trajectories of X1[t] and X2[t] spiral away from the origin.

p1 = ParametricPlot[Evaluate[X1[t]],{t,-3,3}]

p1 = ParametricPlot[Evaluate[X2[t]],{t,-3,3}]

solns = Flatten[Table[X[t],{c1,-2,2},{c2,-2,2}],1];

p = ParametricPlot[Evaluate[solns],{t,-3,3},
    PlotRange->{{-3,3},{-3,3}},AspectRatio->Automatic]

Repeated Eiginvalues (2x2)

Find the general solution of the system X'[t] == A.X[t] where

A = {{ 2, 1},
     {-1, 4}};

Solution

Eigensystem[A]

The  eigenvalue 3 is repeated, and Eigensystem[] finds only one non-trivlial eigenvector, {1,1}. This gives one fundamental solution

X1[t_] = {1,1}E^(3*t)

To find another fundamental solution we need to find a generalized eigenevector. To do this we solve:

Solve[(A - 3*IdentityMatrix[2]).{w1,w2} == {1,1}]

So {-1,0} is one generalized eigenvector giving us a second fundamental solution:

X2[t_] = {1,1} t E^(3*t) + {-1, 0} E^(3*t)

The general solution is:

X[t_] = c1 X1[t] + c2 X2[t]

Let's check:

X'[t]//Simplify

A.X[t]//Simplify

Multiples of the first fundamental solution gives a line:

p1 = ParametricPlot[Evaluate[{X1[t],-X1[t]}], {t,-2,2},
PlotStyle->RGBColor[1,0,0]]

Multiples of the other fundamental solution is a curve tangent to the first line at {0,0}:

p2 = ParametricPlot[Evaluate[{X2[t],-X2[t]}], {t,-2,2},
PlotStyle->RGBColor[0,0,1]]

solns = Flatten[Table[X[t],{c1,-3,3}, {c2,-3,3}],1];

p = ParametricPlot[Evaluate[solns], {t,-2,2},
PlotRange->{{-5,5},{-5,5}}]

Show[p,p1,p2]

Repeated Eiginvalues (3x3)

Find the general solution of the system X'[t] == A.X[t] where

A = {{-4,  2, -1},
     {-1, -4,  2},
     {-1, -1, -1}};

Solution

Eigensystem[A]

The eigenvalue -3 is repeated three times. Only one non-trivial eigenvector, {1,1,1}, is found which gives a fundamental solution:

X1[t_] = {1,1,1}E^(-3t);

We now try to find two generalized eigenvectors. First solve:

Solve[(A-(-3)IdentityMatrix[3]).{w1,w2,w3}=={1,1,1}]

Choose w3 == 0 to get the generalized eigenvector {-1,0,0} and a corresponding fundamental solution:

X2[t_] = {1,1,1}t*E^(-3t) + {-1,0,0}E^(-3t);

To find another, solve:

Solve[(A-(-3)IdentityMatrix[3]).{u1,u2,u3}=={-1,0,0}]

Choose u3 == 0 to get the second generalized eigenvector {1/3,-1/3,0}. A third fundamental solution is then given by

X3[t_] = {1,1,1}(t^2/2)*E^(-3t) + {-1,0,0}t*E^(-3t) +
         {1/3,-1/3,0}E^(-3t);

And the general solution is:

X[t_] = c1 X1[t] + c2 X2[t] + c3 X3[t]//Simplify

Let's check:

X'[t]//Simplify

A.X[t]//Simplify

solns = Flatten[Table[X[t],{c1,-1,1},{c2,-1,1},{c3,-1,1}],2];

p = ParametricPlot3D[Evaluate[solns], {t,-2,2}]

Fundamental Matrix

Find the fundamental matrix Phi[t] such that Phi[0] is the identity matrix for the system
X'[t] == A.X[t]where

A = {{-4, -3, -1},
     {-1,  0,  5},
     { 1,  3, -2}};

Use Phi[t] to find a solution satisfying the given initial condition. X[0]=={1,0,1}.

Solution

Eigensystem[A]

X1[t_] = {-1, -1, 1}E^(-6t);
X2[t_] = { 2, -1, 1}E^(-3t);
X3[t_] = {-1,  2, 1}E^( 3t);

A fundamental matrix has columns given by X1[t], X2[t], and X3[t]. We use Transpose[] to covert these vectors from rows into columns.

Xi[t_] = Transpose[{X1[t], X2[t], X3[t]}];
Xi[t]//MatrixForm

Xi[0]//MatrixForm

To get a fundamental matrix Phi[t] such that Phi[0] is the identity matrix, we multiply Xi[t] on the right by the inverse of Xi[0].

Phi[t_] = Xi[t].Inverse[Xi[0]];
Phi[t]//MatrixForm

Phi[0]//MatrixForm

Let's check that Phi[t] satisfies the equation:

Phi'[t]//Simplify//MatrixForm

A.Phi[t]//Simplify//MatrixForm

The solution satisfying the given initial condition is simply:

X[t_] = Phi[t].{1,0,1}//Simplify

Compare with DSolve[]:

DSolve[
{x1'[t] == -4x1[t]-3x2[t]- x3[t],
 x2'[t] == - x1[t]       +5x3[t],
 x3'[t] ==   x1[t]+3x2[t]-2x3[t],
 x1[0]==1,x2[0]==0,x3[0]==1},
 {x1[t],x2[t],x3[t]},t]//Simplify

Variation of Parameters

Use the method of variation of parameters to find the general solution of the non-homogeneous system
X'[t] == A.X[t] + G[t], X[0] == {1,0,1},
where

A = {{-4, -3, -1},
     {-1,  0,  5},
     { 1,  3, -2}};
G[t_] = {1,t,t^2};

Solution

A solution to the associated homogeneous problem,  X'[t] == A.X[t], was found in in the previous problem:

X1[t_] = {-1, -1, 1}E^(-6t);
X2[t_] = { 2, -1, 1}E^(-3t);
X3[t_] = {-1,  2, 1}E^( 3t);
Xi[t_] = Transpose[{X1[t],X2[t],X3[t]}];
IXi[t_] = Inverse[Xi[t]]//Simplify;
IXi[t]//MatrixForm

A particular solution to the non-homogeneous problem is given by Y[t_] = Xi[t].U[t] where U[t] satisfies U'[t] == IXi[t].G[t].

IXi[t].G[t]

U[t_] = Integrate[IXi[t].G[t], t]//Simplify

Y[t_] = Xi[t].U[t]//Simplify

The general solution to the given system is then

X[t_] = c1 X1[t] + c2 X2[t] + c3 X3[t] + Y[t]

Let's check:

X'[t]//Simplify

A.X[t] + G[t]//Simplify

We now solve for the constants c1, c2, and c3 that satisfy the given initial condition.

cs = Solve[X[0] == {1,0,1}]

The solution to the initial value problem is thus

X[t_] = X[t] /. First[cs] //Simplify

The solution given by DSolve[] is:

DSolve[
{x1'[t] == -4x1[t]-3x2[t]- x3[t] + 1,
 x2'[t] == - x1[t]       +5x3[t] + t,
 x3'[t] ==   x1[t]+3x2[t]-2x3[t] + t^2,
 x1[0]==1,x2[0]==0,x3[0]==1},
 {x1[t],x2[t],x3[t]},t]//Simplify