Part IV. Linear Differential Equations

Section 1. Linear Oscillators

2.* Consider now the following damped system. Obtain the solution.

y'' + y' + 4y = 0  , y(0) = 2 , y'(0) = -3

From the form of the solution decide if the system is underdamped, critically damped, or overdamped.    

In[93]:=

DE = y''[t] + y'[t] + 4y[t] == 0
soln = DSolve[ {DE, y[0]==2, y'[0]==-3}, y, t]

Out[93]=

4 y[t] + y^′[t] + y^′′[t] 0

Out[94]=

{{yFunction[{t}, 2/15 ^(-t/2) (15 Cos[(15^(1/2) t)/2] - 2 15^(1/2) Sin[(15^(1/2) t)/2])]}}

The system is underdamped.

a.   What is the pseudo-period of the oscillations?   

The pseudo-period is

In[95]:=

N[2Pi/(Sqrt[15]/2)]*time_units

Out[95]=

3.24462 time_units

b.   What is the time constant?   

The time constant is 2 time units.

c.   Based upon your answer to part b estimate the time interval required for the oscillations to disappear  from view.   

In 10 time units the oscillations will decrease by a factor of e^5 which is about 150. They will disappear.

d.   Plot the solution curve over the interval you named in part c.   

In[96]:=

Plot[ y[t]/.soln, {t,0,10}, PlotRange->{-1.5,2}]

[Graphics:../HTMLFiles/Solutions_93.gif]

e.   Add to the curve in part d the curves defined by A e^(-t/2)and   -A e^(-t/2)where  A = (4 + 16/15)^(1/2)  .   Make them blue. What is the significance of these curves? Where did the formula for A come from?   

In[97]:=

A = Sqrt[4 + 16/15];
Plot[ {y[t]/.soln,A Exp[-t/2],-A Exp[-t/2]}, {t,0,10}, PlotRange->{-2.2,2.2}, PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1],RGBColor[0,0,1]}]

[Graphics:../HTMLFiles/Solutions_97.gif]

The blue curves envelope the solution curve. The constant A is derived from the standard formula used to convert the sum of a sine and a cosine into a sinusoid.

f.    Convert the solution into the function g. Use g to plot the phase plane trajectory.    

In[101]:=

g = soln[[1,1,2]];
g[t]

Out[102]=

2/15 ^(-t/2) (15 Cos[(15^(1/2) t)/2] - 2 15^(1/2) Sin[(15^(1/2) t)/2])

In[113]:=

traj = ParametricPlot[ {g[t],g'[t]}, {t,0,10}, PlotRange->{{-2,2},{-4,2}}, AspectRatio->1/1]

[Graphics:../HTMLFiles/Solutions_99.gif]

g.    Add to the trajectory the points corresponding to t = 0, 0.25, 0.5, 0.75, ... , 2.0.

In[117]:=

Show[ ListPlot[ Table[{g[t],g'[t]},{t,0,2,0.25}],
                PlotStyle->PointSize[0.02]],
      traj, PlotRange->{{-2,2},{-4,2}}, AspectRatio->1/1]

[Graphics:../HTMLFiles/Solutions_100.gif]

Section 2. State Space

2.* Damp the system slightly by changing the IVP to the following

y'' + 0.1y' + 4y = cos(1.8 t)  , y(0) = 2 , y'(0) = -3

Obtain the solution, call it soln and convert it into a function h.     

In[165]:=

DE = y''[t] + 0.1y'[t] + 4y[t] == Cos[1.8t]
soln = DSolve[ {DE,y[0]==2,y'[0]==-3},y,t];
h = soln[[1,1,2]];
N[h[t]]//FullSimplify

Out[165]=

4 y[t] + 0.1 y^′[t] + y^′′[t] Cos[1.8 t]

Out[168]=

(1.2459 - 5.55112*10^-17 ) Cos[1.8 t] + (0.295082 + 2.77556*10^-17  ... 5112*10^-17 ) Cos[1.99937 t] - (1.74727 - 1.11022*10^-16 ) Sin[1.99937 t])

a.    Based upon the solution formula determine the time constant for the beats. That is, how long will it take (approximately) for the beats to disappear from the solution curve as it settles down to its  steady-state mode?   

The time constant for the transient solution (when the beats will appear) is 1/0.05 = 20 time units.

b.    Plot the solution curve to verify your answer to part a.    

The beats should disappear in about 100 time units.

In[145]:=

Plot[ h[t], {t,0,100}, AspectRatio->1/4]

[Graphics:../HTMLFiles/Solutions_104.gif]

c.    Use h to obtain the phase plane and state space trajectories for this system.

In[150]:=

ParametricPlot[ {h[t],h'[t]}, {t,0,100}, PlotLabel->"Phase Plane Trajectory", AspectRatio->3/2]

[Graphics:../HTMLFiles/Solutions_105.gif]

In[177]:=

ParametricPlot3D[ {h[t],h'[t],t}, {t,0,100}, PlotLabel->"State Space Trajectory", AspectRatio→1/1, PlotPoints→400, ViewPoint->{2.664, -1.611, 10.327} ]

[Graphics:../HTMLFiles/Solutions_106.gif]

4.* Obtain a numeric solution for the IVP in Exercise 2. Sketch the time series for position and the time series for velocity. Then plot the phase plane trajectory and the  state space trajectory.

In[181]:=

soln = NDSolve[ {DE,y[0]==2,y'[0]==-3},y,{t,0,100}];
h = soln[[1,1,2]];

In[185]:=

Plot[ h[t], {t,0,100}, AspectRatio->1/3, PlotLabel->"Time Series, Position"]

[Graphics:../HTMLFiles/Solutions_107.gif]

In[186]:=

Plot[ h'[t], {t,0,100}, AspectRatio->1/3, PlotLabel->"Time Series, Velocity"]

[Graphics:../HTMLFiles/Solutions_108.gif]

In[187]:=

ParametricPlot[ {h[t],h'[t]}, {t,0,100}, PlotLabel->"Phase Plane Trajectory", AspectRatio->3/2]

[Graphics:../HTMLFiles/Solutions_109.gif]

In[188]:=

ParametricPlot3D[ {h[t],h'[t],t}, {t,0,100}, PlotLabel->"State Space Trajectory", AspectRatio→1/1, PlotPoints→400, ViewPoint->{2.664, -1.611, 10.327} ]

[Graphics:../HTMLFiles/Solutions_110.gif]

Section 3. Two Dimensional Systems

In[1]:=

<<Graphics`PlotField`

4.* The model

x' = x (1 - y - x/a)
y' = y (1 - x - y/b)

is used to study competing species. See Ledder, Section 5.5, Exercise 11. Use PlotVectorField and DSolve to do the  following.   

a.   Draw the direction field when 1/a = 1.9 and 1/b = 1.5. Use the window 0 ≤  x ≤ 1, 0 ≤ y ≤ 1.   

In[9]:=

a = 1/1.9; b = 1/1.5;
u = x(1-y-x/a); v = y(1-x-y/b);
vf= PlotVectorField[ {u,v}/Sqrt[u^2+v^2], {x,0,1}, {y,0,1}, Axes->True]

                                     1 Power :: infy : Infinite expression  - encountered. More…                                      0

∞ :: indet : Indeterminate expression 0 ComplexInfinity encountered. More…

∞ :: indet : Indeterminate expression 0 ComplexInfinity encountered. More…

[Graphics:../HTMLFiles/Solutions_114.gif]

Out[11]=

⁃Graphics⁃

b.   Find the stationary point and add the nullclines to the direction field. Discuss the solutions based upon the picture you see.   

In[5]:=

statpt = Solve[ {u==0,v==0}, {x,y}]

Out[5]=

{{x0., y0.}, {x0., y0.666667}, {x0.27027, y0.486486}, {x0.526316, y0.}}

In[6]:=

SP = {x,y}/.statpt[[3]]

Out[6]=

{0.27027, 0.486486}

In[43]:=

<<Graphics`ImplicitPlot`
nullclines = ImplicitPlot[{u==0,v==0}, {x,0,1}, {y,0,1}];
stationarypoint = ListPlot[ {SP}, PlotStyle->PointSize[0.025]];
flow = Show[ nullclines, stationarypoint, vf]

[Graphics:../HTMLFiles/Solutions_118.gif]

Based on this picture, all solution trajectories in the first quadrant will move towards the stationary point.

c.   Add solution curves corresponding to the following set of initial conditions (a circle of points around  the stationary point)   x(0) = 0.3 + 0.2cos(k π/6) , y(0) = 0.5 + 0.2sin(k π/6) , k = 1, 2, ... , 12

In[38]:=

solns = Table[
  NDSolve[ {x'[t]==x[t](1-y[t]-x[t]/a),
            y'[t]==y[t](1-x[t]-y[t]/b),
                 x[0]==0.3+0.2Cos[k*Pi/6],y[0]==0.5+0.2Sin[k*Pi/6]},
                 {x[t],y[t]}, {t,-10,10}], {k,11}];

In[42]:=

Show[ flow,
      Table[ ParametricPlot[ {x[t],y[t]}/.solns[[k]], {t,-10,10}, PlotStyle->{RGBColor[1,0,0],Thickness[0.01]} ], {k,11}],
      PlotRange->{{0,1},{0,1}}]

[Graphics:../HTMLFiles/Solutions_120.gif]

Section 4. Matrix Methods

2.* Define the matrix B having the columns {0,1,-1}, {1,1,0}, {-1,0,1}. Hint.

B = Transpose[ { {0,1,-1},  {1,1,0}, {-1,0,1} } ]

In[48]:=

B = Transpose[ {{0,1,-1},{1,1,0},{-1,0,1}}];
B//MatrixForm

Out[49]//MatrixForm=

( 0    1    -1 )            1    1    0            -1   0    1

a.   Find the characteristic polynomial of B and factor it.   

In[50]:=

CharacteristicPolynomial[B,t]
Factor[%]

Out[50]=

-2 + t + 2 t^2 - t^3

Out[51]=

-(-2 + t) (-1 + t) (1 + t)

b.   Find B's eigenvalues with Eigenvalues[B]. Name them lambda.

In[52]:=

lambda = Eigenvalues[B]

Out[52]=

{2, -1, 1}

c.   Find B's eigenvectors using Eigenvectors[B]. Name them V and define P to be the transpose of V.   

In[58]:=

V = Eigenvectors[B]
P = Transpose[V]
P//MatrixForm

Out[58]=

{{-1, -1, 1}, {2, -1, 1}, {0, 1, 1}}

Out[59]=

{{-1, 2, 0}, {-1, -1, 1}, {1, 1, 1}}

Out[60]//MatrixForm=

( -1   2    0  )            -1   -1   1            1    1    1

d.   Calculate  P^(-1) B P

In[61]:=

Inverse[P].B.P//MatrixForm

Out[61]//MatrixForm=

( 2    0    0  )            0    -1   0            0    0    1

3.* Continuing 2. Obtain the solution to v' = Bv satisfying v(0) = {1,2,3}}.   

a.    Do if first using DSolve applied to the system of linear differential equations defined by v' = Bv with  the appropriate initial conditions.   

In[68]:=

v[t_] := {x[t],y[t],z[t]}
DEsys = Table[ v'[t][[k]] == (B.v[t])[[k]], {k,3}]

Out[69]=

{x^′[t] y[t] - z[t], y^′[t] x[t] + y[t], z^′[t]  -x[t] + z[t]}

In[72]:=

soln1 = DSolve[ Join[DEsys,{x[0]==1,y[0]==2,z[0]==3}], {x,y,z}, t]

Out[72]=

{{xFunction[{t}, ^(-t)], yFunction[{t}, 1/2 ^(-t) (-1 + 5 ^(2 t))], zFunction[{t}, 1/2 ^(-t) (1 + 5 ^(2 t))]}}

Check:

In[75]:=

DEsys/.soln1//FullSimplify

Out[75]=

{{True, True, True}}

b.   Do it second by making a fundamental matrix solution X(t) defined as the matrix with the eigenvector  solutions in the columns then computing    

v(t) = X(t) X(0)^(-1) v(0)

In[79]:=

X = Transpose[Table[ Exp[lambda[[k]]t]*V[[k]], {k,3}]]

Out[79]=

{{-^(2 t), 2 ^(-t), 0}, {-^(2 t), -^(-t), ^t}, {^(2 t), ^(-t), ^t}}

In[81]:=

vsoln = X.Inverse[X/.t->0].{1,2,3};
vsoln//FullSimplify

Out[82]=

{^(-t), 2 Cosh[t] + 3 Sinh[t], 3 Cosh[t] + 2 Sinh[t]}

That is interesting output. Let's check that the solution displayed above is the same as the solution we got in part a.

In[83]:=

{x[t],y[t],z[t]}/.soln1//FullSimplify

Out[83]=

{{^(-t), 2 Cosh[t] + 3 Sinh[t], 3 Cosh[t] + 2 Sinh[t]}}

Indeed, it is.

c.   Do it third by using the Matrix Exponential, Mexp. Once you have it, the solution is  

v(t) = Mexp v(0)  

In[92]:=

Mexp = MatrixExp[B*t];
Mexp//MatrixForm

Out[93]//MatrixForm=

( 1         -t              3 t                       1         -t               3 t   ...                                               6                                                  6

In[94]:=

vsoln2 = Mexp.{1,2,3};
vsoln2//FullSimplify

Out[95]=

{^(-t), 2 Cosh[t] + 3 Sinh[t], 3 Cosh[t] + 2 Sinh[t]}

Section 5. The Laplace Transform

5.* Consider the following initial value problem

y'' + y' + 9 y = UnitStep(t - 2 π ) - UnitStep(t - 5 π ) + 2 DiracDelta(t - 9 π )  ,  y(0) = 0 , y'(0) = 1

a. Obtain the solution using DSolve. Plot it for 0 ≤ t ≤ 40  and explain the behavior of the solution curve.  

In[1]:=

DE = y''[t] + y'[t] + 9y[t] == UnitStep[t-2Pi] - UnitStep[t-5Pi] + 2DiracDelta[t-9Pi]
soln = DSolve[ {DE,y[0]==0,y'[0]==1}, y, t];

Out[1]=

9 y[t] + y^′[t] + y^′′[t] 2 DiracDelta[-9 π + t] - UnitStep[-5 π + t] + UnitStep[-2 π + t]

In[3]:=

Plot[ y[t]/.soln, {t,0,40}, PlotRange->{-0.4,0.6}, AspectRatio->1/3]

[Graphics:../HTMLFiles/Solutions_140.gif]

The solution curve shows the response to the steady force from t = 6 to t = 16 (roughly). See the plot of  the driver below. Then there is no force until the kick provided by the Dirac delta at about t = 27. Note  that at t = 27 the system reacts as it did at t = 0 except it goes twice as far from equilibrium.

In[4]:=

Plot[ DE[[2]], {t,0,40}, PlotRange->{-0.1,1.4}, PlotLabel->"Unit Step Driver", AspectRatio->1/3]

[Graphics:../HTMLFiles/Solutions_141.gif]

b.   Obtain the solution using the method of Laplace transforms. How does the solution formula compare to the  formula obtained in part a?

The Laplace transform method:
Transform the DE to TDE (in the s-domain)

In[5]:=

TDE = LaplaceTransform[DE, t, s]

Out[5]=

9 LaplaceTransform[y[t], t, s] + s LaplaceTransform[y[t], t, s] + s^2 LaplaceTransform[y[t], t ... ^′[0] 2 ^(-9 π s) - ^(-5 π s)/s + ^(-2 π s)/s

Substitute the initial values and introduce Y for the transform of the solution.

In[6]:=

TDE/.{y[0]->0,y'[0]->1,LaplaceTransform[y[t],t,s]->Y}

Out[6]=

-1 + 9 Y + s Y + s^2 Y2 ^(-9 π s) - ^(-5 π s)/s + ^(-2 π s)/s

Solve for Y

In[7]:=

Y = Y/.Simplify[Solve[%,{Y}]][[1]]

Out[7]=

(^(-9 π s) (-^(4 π s) + ^(7 π s) + 2 s + ^(9 π s) s))/(s (9 + s + s^2))

Invert to get the Laplace transform solution.

In[8]:=

Lsoln = InverseLaplaceTransform[Y,s,t]

Out[8]=

(2 ^(-t/2) Sin[(35^(1/2) t)/2])/35^(1/2) + (4 ^(1/2 (9 π - t)) Sin[1/2 35 ... 2 35^(1/2) (-2 π + t)] + 35^(1/2) Sin[1/2 35^(1/2) (-2 π + t)])) UnitStep[-2 π + t]

This is pretty complicated. The solution using DSolve is even worse.

In[11]:=

y[t]/.soln

Out[11]=

{-1/315 ^(-t/2) (-18 35^(1/2) Sin[(35^(1/2) t)/2] + 36 35^(1/2) ^(9 π/2)  ...  t)/2] UnitStep[-2 π + t] - 35 ^(t/2) Sin[(35^(1/2) t)/2]^2 UnitStep[-2 π + t])}

6.* Use one of the solution formulas in 5 a to make a function, g. Use g to plot the phase plane trajectory  and as well as the state space trajectory for the system.

In[13]:=

g[t_] := Evaluate[Lsoln]
g[t]

Out[14]=

(2 ^(-t/2) Sin[(35^(1/2) t)/2])/35^(1/2) + (4 ^(1/2 (9 π - t)) Sin[1/2 35 ... 2 35^(1/2) (-2 π + t)] + 35^(1/2) Sin[1/2 35^(1/2) (-2 π + t)])) UnitStep[-2 π + t]

In[20]:=

ParametricPlot[ {g[t],g'[t]}, {t,0,40}, PlotRange->{{-0.4,0.6},{-1.5,2.2}}, AspectRatio->3/2]

[Graphics:../HTMLFiles/Solutions_148.gif]

In[40]:=

ParametricPlot3D[ {g[t],g'[t],t}, {t,0,40}, AspectRatio→3/2, ViewPoint->{3, -15, 50}, PlotPoints→300, PlotRange->{{-0.4,0.6},{-1.5,2.2},{0,40}}]

[Graphics:../HTMLFiles/Solutions_149.gif]

7.* Continuing 6. Plot an animation of the phase space trajectory.

In[41]:=

<<Graphics`Animation`

In[42]:=

Animate[ ParametricPlot[ {g[t],g'[t]}, {t,0,T}, PlotRange->{{-0.4,0.6},{-1.5,2.2}}, AspectRatio->3/2], {T,1,40}]

[Graphics:../HTMLFiles/Solutions_174.gif]


Created by Mathematica  (December 8, 2004)