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]=
Out[94]=
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]=
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 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}]
e. Add to the curve in part d the curves defined by and
where A =
. 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]}]
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]=
In[113]:=
traj = ParametricPlot[ {g[t],g'[t]}, {t,0,10}, PlotRange->{{-2,2},{-4,2}}, AspectRatio->1/1]
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]
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]=
Out[168]=
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 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]
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]
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} ]
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"]
In[186]:=
Plot[ h'[t], {t,0,100}, AspectRatio->1/3, PlotLabel->"Time Series, Velocity"]
In[187]:=
ParametricPlot[ {h[t],h'[t]}, {t,0,100}, PlotLabel->"Phase Plane Trajectory", AspectRatio->3/2]
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} ]
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]
Out[11]=
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]=
In[6]:=
SP = {x,y}/.statpt[[3]]
Out[6]=
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]
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)
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}}]
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=
a. Find the characteristic polynomial of B and factor it.
In[50]:=
CharacteristicPolynomial[B,t]
Factor[%]
Out[50]=
Out[51]=
b. Find B's eigenvalues with Eigenvalues[B]. Name them lambda.
In[52]:=
lambda = Eigenvalues[B]
Out[52]=
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]=
Out[59]=
Out[60]//MatrixForm=
d. Calculate
In[61]:=
Inverse[P].B.P//MatrixForm
Out[61]//MatrixForm=
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]=
In[72]:=
soln1 = DSolve[ Join[DEsys,{x[0]==1,y[0]==2,z[0]==3}], {x,y,z}, t]
Out[72]=
Check:
In[75]:=
DEsys/.soln1//FullSimplify
Out[75]=
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)
In[79]:=
X = Transpose[Table[ Exp[lambda[[k]]t]*V[[k]], {k,3}]]
Out[79]=
In[81]:=
vsoln = X.Inverse[X/.t->0].{1,2,3};
vsoln//FullSimplify
Out[82]=
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]=
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=
In[94]:=
vsoln2 = Mexp.{1,2,3};
vsoln2//FullSimplify
Out[95]=
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]=
In[3]:=
Plot[ y[t]/.soln, {t,0,40}, PlotRange->{-0.4,0.6}, AspectRatio->1/3]
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]
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]=
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]=
Solve for Y
In[7]:=
Y = Y/.Simplify[Solve[%,{Y}]][[1]]
Out[7]=
Invert to get the Laplace transform solution.
In[8]:=
Lsoln = InverseLaplaceTransform[Y,s,t]
Out[8]=
This is pretty complicated. The solution using DSolve is even worse.
In[11]:=
y[t]/.soln
Out[11]=
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]=
In[20]:=
ParametricPlot[ {g[t],g'[t]}, {t,0,40}, PlotRange->{{-0.4,0.6},{-1.5,2.2}}, AspectRatio->3/2]
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}}]
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}]
Created by Mathematica (December 8, 2004)