Math 325: Differential Equations	Solutions to Assignment 3		Fall 1997

Initialization

<<Graphics`PlotField`

Problem 1

Using eigenvalues and eigenvectors,  express the general solution of the system X'[t] == A.X[t] in terms of real-valued functions for each of the matrices A below. Check that your solutions are correct by comparing X'[t] and A.X[t]. Superimpose several plots for each using ParametricPlot3D[] or ParametricPlot3D[].

a)  A = {{ 3, 2},
       {-1, 1}};

b)  A = {{1,-1, 4},
       {3, 2,-1},
       {2, 1,-1}};
       
c)  A = {{-14,18, 8},
       { -6, 8, 2},
       { -5, 7, 3}};

Solution

a)

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

Eigensystem[A]

A complex-valued solution is:

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

The real and imaginary parts of Z[t] give two indepenedent real-valued solutions:

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 general solution is

X'[t]//Simplify

A.X[t]//Simplify

Now we plot several solutions.

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

The output of Table[] is a doubly-nested list; Flatten[..,1] removes the outer layer leaving a list of vector functions. For example, the first on the list would be X[t] when both constants are -1:

solns[[1]]

ParametricPlot[Evaluate[solns], {t,-1,1}]

The solutions are stable spirals that flow into the origin.

b)

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

Eigensystem[A]

Fundamental solutions are:

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

The general solution is

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

Check:

X'[t]//Simplify

A.X[t]//Simplify

Plots:

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

In this case, Table[] creates a triply-nested list; Flatten[..,2] removes the outer two layers.

ParametricPlot3D[Evaluate[solns], {t,-1,1}]

c)

A = {{-14,18, 8},
     { -6, 8, 2},
     { -5, 7, 3}};

Eigensystem[A]

We get a real-valued solution and a complex-valued solution:

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

The real and imaginary parts of Z[t] give two indepenedent real-valued solutions:

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

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

The general solution is

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

Check:

X'[t]//Simplify

A.X[t]//Simplify

Plots:

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

ParametricPlot3D[Evaluate[solns], {t,-Pi/2,Pi/2},
ViewPoint->{1.323, -2.716, -1.525},
PlotRange->{{-5,5},{-5,5},{-2,2}}]

Problem 2

For each of the matrices A below, find the fundamental matrix Phi[t] such that Phi[0] is the identity matrix for the system X'[t] == A.X[t]. Verify that Phi'[t]==A.Phi[t]. Then use Phi[t] to find a solution satisfying the given initial condition.

a)  A = {{-7,-8, 2},
        { 9,11,-3},
        { 9, 8, 0}};
  X[0]=={3,-1,1};

b)  A = {{-4,-2,-1},
        { 2, 6,-3},
        { 2, 4,-1}};
  X[0]=={1,1,-2};
 
c)  A = {{  5, 3,-1},
        {-16,-6,-1},
        { -4,-3, 2}};
  X[0]=={1,1,1};

Solution

a)

A = {{-7,-8, 2},
     { 9,11,-3},
     { 9, 8, 0}};

Eigensystem[A]

Fundamental solutions:

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

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

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

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

Phi[t_] = Psi[t].Inverse[Psi[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].{3,-1,1}

b)

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

Eigensystem[A]

We get two fundamental solutions from this:

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

The eigenvalue 2 is repeated and Eigensystem[] finds only one non-trivial eigenvector, {-1,2,2}. A generalized eigenvector is found by solving

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

For example, {w1,w2,w3} = {-1,5/2,2}, from which we can derive a third independent solution:

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

The general solution is:

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

Check:

X'[t]//Expand

A.X[t]//Expand

Fundamental matrix:

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

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

The solution satisfying the given initial condition is:

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

c)

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

Eigensystem[A]

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

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

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

General solution:

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

Check:

X'[t]//Simplify

A.X[t]//Simplify

Fundamental matrix:

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

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

The solution satisfying the given initial condition is:

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

Problem 3

Use the method of variation of parameters to find the general solution of the system X'[t] == A.X[t] + G[t] for each of the matrices A and vector functions G[t] below. Check that your solution is correct by comparing X'[t] and A.X[t]+G[t]. Then find the solution that satisfies the given initial condition.

a) (see Problem 2a)
   A = {{-7,-8, 2},
        { 9,11,-3},
        { 9, 8, 0}};
  G[t_] = {t^2,t,1};
  X[0]=={3,-1,1};

b) (see Problem 2b)
   A = {{-4,-2,-1},
        { 2, 6,-3},
        { 2, 4,-1}};
  G[t_] = {1,Cos[t],Sin[t]};
  X[0]=={1,1,-2};
 
c) (see Problem 2c)
   A = {{  5, 3,-1},
        {-16,-6,-1},
        { -4,-3, 2}};
  G[t_] = {E^t,1,E^(-t)};
  X[0]=={1,1,1};

Solution

a)

A = {{-7,-8, 2},
     { 9,11,-3},
     { 9, 8, 0}};
G[t_] = {t^2,t,1};

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

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

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

IPsi[t].G[t]

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

Y[t_] = Psi[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]

Check:

X'[t]//Simplify

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

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

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

The solution to the initial value problem is thus

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

b)

A = {{-4,-2,-1},
     { 2, 6,-3},
     { 2, 4,-1}};
G[t_] = {1,Cos[t],Sin[t]};

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

X1[t_] = {-3,1,1}E^(-3t);
X2[t_] = {-1,2,2}E^(2t);
X3[t_] = {-1,2,2} t*E^(2t) + {-1,5/2,2} E^(2t);
Psi[t_] = Transpose[{X1[t], X2[t], X3[t]}];
IPsi[t_] = Inverse[Psi[t]]//Simplify;
IPsi[t]//MatrixForm

A particular solution is Y[t_] = Psi[t].U[t] where U'[t] == IPsi[t].G[t].

IPsi[t].G[t]//Simplify

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

Y[t_] = Psi[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]

Check:

X'[t]//Simplify

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

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

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

The solution to the initial value problem is thus

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

c)

A = {{  5, 3,-1},
     {-16,-6,-1},
     { -4,-3, 2}};
G[t_] = {E^t,1,E^(-t)};

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

X1[t_] = {-1,2,2}E^t;
X2[t_] = {-Cos[3 t],2 Cos[3 t]+Sin[3 t],Cos[3 t]};
X3[t_] = {-Sin[3 t],-Cos[3 t]+2 Sin[3 t],Sin[3 t]};
Psi[t_] = Transpose[{X1[t], X2[t], X3[t]}];
IPsi[t_] = Inverse[Psi[t]]//Simplify;
IPsi[t]//MatrixForm

IPsi[t].G[t]//Simplify

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

Y[t_] = Psi[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]

Check:

X'[t]//Simplify

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

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

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

The solution to the initial value problem is thus

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

Problem 4

Consider the autonomous system:
    x'[t] == 1 - x[t]y[t]
    y'[t] == x[t] - y[t]^3
a)  Find all critical points of this system. At each critical point, calculate the corresponding linear system and find the eigenvalues of the coefficient matrix. Then identify the type and stability of the critical point.
b)  Plot the vector field using PlotVectorField[] with the options ScaleFunction->(Log[50#+1]&) (this makes the arrows bigger near the critical points) and Frame->True. Use the vector field plot to choose a representative sample of several initial points {x0,y0}. Use NDSolve[] to obtain numerical solutions through these points and plot the corresponding trajectories on a single graph. Indicate what happens to the trajectories for increasing t.

Solution

F[x_, y_] = 1 - x y;
G[x_, y_] = x - y^3;

a)  Find the critical points:

cp = Solve[{F[x,y] == 0, G[x,y] == 0}]

The critical points are thus {-1,-1}, and {1,1}. To analyze them we need to examine the matrix of partial derivative evaluated at each.

A = {{D[F[x,y],x], D[F[x,y],y]},
     {D[G[x,y],x], D[G[x,y],y]}};
A//MatrixForm

At the critical point {-1,-1}:

A1 = A /. cp[[1]];
A1//MatrixForm

Eigensystem[A1]

N[%]

Since the eigenvalues have opposite signs, the point {-1,-1} is unstable saddle point.
At the point {1,1}:

A2 = A /. cp[[4]];
A2//MatrixForm

Eigensystem[A2]

Since there is only one independent eigenvector, the point {1,1} is a node (proper or improper) or a spiral point. It is asymtotically stable since both eigenvalues are negative.

b)

?PlotVectorField

v = PlotVectorField[{F[x,y],G[x,y]},{x,-2,2},{y,-2,2},
ScaleFunction->(Log[50#+1]&), Frame->True]

From the vector field plot, it appears that initial points chosen around the edge of the above rectangle will produce a representative set of trajectories. Both First[] and Flatten[] are present below only to get rid of extra braces. Evaluate[] keeps Mathematica from complaining about not being able to compile the functions.

solns1 = Flatten[Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == 1 - x[t]y[t],
         y'[t] == x[t] - y[t]^3,
         x[0] == i, y[0] == j}, {x[t],y[t]},{t,0,10}]],
        {i,-2,2,4}, {j,-2,2,.5}],1];

solns2 = Flatten[Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == 1 - x[t]y[t],
         y'[t] == x[t] - y[t]^3,
         x[0] == i, y[0] == j}, {x[t],y[t]},{t,0,10}]],
        {i,-2,2,.5}, {j,-2,2,4}],1];

p = ParametricPlot[Evaluate[Union[solns1,solns2]],{t,0,10},
PlotRange->{{-4,4},{-4,4}}]

It appears that the trajectories near and to the left of the saddle point  at {-1,-1} go off to infinity asymptotically along the same curve towards the left. Trajectories to the right of the saddle point seem to spiral into the critical point at {1,1}.

Show[p,v,PlotRange->{{-2,2},{-2,2}}]

Problem 5

In this problem we consider the equations describing the motion of a pendulum, both damped and undamped. For each situation below convert the equation into a system of equations in x[t] = theta[t] and y[t] = theta'[t]. Plot the corresponding vector field. Use NDSolve[] to solve this system and plot on a single graph the trajectories corresponding to the given initial velocities b. Based on these plots, describe what the pendulum is doing in each case. Experiment to see if you can find a value for the initial velocity b such that the trajectory approaches the unstable equilibrium point {Pi,0} corresponding to the pendulum being stationary in an upright position.

a) The undamped pendulum:

    theta''[t] + Sin[theta] == 0,
    theta[0] == 0, theta'[0] == b

The initial velocities are b = 0.5 to 3.0 in steps of 0.5.

b)  The damped pendulum:

    theta''[t] + 0.5 theta'[t] + Sin[theta] == 0,
    theta[0] == 0, theta'[0] == b

The initial velocities are b = 0.5 to 6.0 in steps of 0.5.

Solution

a)  Since x'[t] == theta'[t] == y[t], and y'[t]==theta''[t]==-Sin[x[t]], the system is:

    x'[t] == y[t],       x[0] == 0
    y'[t] == -Sin[x[t]], y[0] == b

PlotVectorField[{y,-Sin[x]},{x,-2,16},{y,-4,4},
    ScaleFunction->(Log[50#+1]&),Frame->True]

solns = Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],
         y'[t] == -Sin[x[t]],
         x[0] == 0, y[0] == b}, {x[t],y[t]}, {t,0,20}]],
        {b, 0.5, 3, 0.5}];

ParametricPlot[Evaluate[solns],{t,0,20},PlotRange->{{-2,16},{-2,4}},
AspectRatio->1/3]

The first three initial velocities start the pendulum swinging back and forth so that the same x (theta) and y (theta') coordinates are repeated over and over. This translates into the three closed trajectories in the xy-plane. The initial velocities 2, 2.5, and 3.0 are enough to make the pendulum go all the way around, so that the angle (x==theta) is ever increasing. There is an initial velocity which lies on the boundary of these two types of behaviour: namely an initial velocity just enough to swing the pendulum from pointing down (x==0) to pointing up (x==Pi). The trajectory would approach the unstable equilibrium point {Pi,0}. We know that an initial velocity of 2 is too fast. Let's experiment with a smaller values to see if we can `shrink' the trajectory so that it stops at the point {Pi,0}.

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
    {x'[t] == y[t],y'[t] == -Sin[x[t]],x[0] == 0, y[0] == 1.99},
    {x[t],y[t]},{t,0,20}]]],{t,0,20},AspectRatio->Automatic]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
    {x'[t] == y[t],y'[t] == -Sin[x[t]],x[0] == 0, y[0] == 1.9999},
    {x[t],y[t]},{t,0,20}]]],{t,0,20},AspectRatio->Automatic]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
   {x'[t] == y[t],y'[t] == -Sin[x[t]],x[0] == 0, y[0] == 1.9999995},
   {x[t],y[t]},{t,0,20}]]],{t,0,20},AspectRatio->Automatic]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
 {x'[t] == y[t],y'[t] == -Sin[x[t]],x[0] == 0, y[0] == 1.999999518},
 {x[t],y[t]},{t,0,20}]]],{t,0,20},AspectRatio->Automatic]

It looks like an initial velocity of 1.999999518 is close to the critical value. Of course, numerical methods are being used, so we should be skeptical of any such numbers until we make a closer analysis of the error.

b) The system for the damped system is deriveed as in part a:

    x'[t] == y[t],                 x[0] == 0
    y'[t] == -Sin[x[t]] - 0.5 y[t], y[0] == b

PlotVectorField[{y,-Sin[x]-0.5y},{x,-2,10},{y,-2,6},
    ScaleFunction->(Log[50#+1]&), Frame->True]

solns = Table[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],
         y'[t] == -Sin[x[t]] - 0.5 y[t],
         x[0] == 0, y[0] == b}, {x[t],y[t]}, {t,0,20}]],
        {b, 0.0, 6.0, 0.5}];

ParametricPlot[Evaluate[solns],{t,0,20},PlotRange->{{-1,10},{-2,7}}]

Here the critical velocity appears to be around 3.1. For initial velocities below this, the pendulum swings back and forth with ever decreasing angles x (theta) and velocity y (theta') coordinates. For initial velocities above 3.1, the pendulum swings around once, bringing all the angles to a range between Pi and 3Pi. The damping causes the pendulum to slow down and again the velocity decreases to 0 while the the angles approach 2Pi. Let's plot some solutions near the critical velocity.

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -Sin[x[t]]-0.5y[t],
         x[0] == 0, y[0] == 3.1},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},PlotRange->{{-1,8},{-2,4}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -Sin[x[t]]-0.5y[t],
         x[0] == 0, y[0] == 3.08},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},PlotRange->{{-1,8},{-2,4}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -Sin[x[t]]-0.5y[t],
         x[0] == 0, y[0] == 3.0885},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},PlotRange->{{-1,8},{-2,4}}]

ParametricPlot[Evaluate[{x[t],y[t]} /. First[NDSolve[
        {x'[t] == y[t],y'[t] == -Sin[x[t]]-0.5y[t],
         x[0] == 0, y[0] == 3.0885098},
        {x[t],y[t]},{t,0,20}]]],{t,0,20},PlotRange->{{-1,8},{-2,4}}]

It looks like an initial velocity of 3.0885098 is close to the critical value.