Part IV. Linear Differential Equations

Section 1. Linear Oscillators

The harmonic oscillator is the fundamental model for the analysis of oscillating systems. Phase plane  trajectories are constructed. PlotVectorField draws direction fields.

1. Obtain the solution to the following initial value problem. Call it soln.  

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

From the form of the solution decide if the system is undamped, underdamped, critically damped, or  overdamped. What is the period of the oscillations?

2, (Continuing 1) Plot the solution to the IVP in Exercise 1. Then create a plot showing the cosine term, the  sine term and their sum (the solution curve). Make the solution red, the cosine blue, and the sine green.  Hint. The cosine term is the solution using the initial conditions y(0) = 2, y'(0) = 0.

a.  What IVP does the sine term solve?   
b.  Determine the amplitude of the oscillations by solving y'(t) = 0 and substituting the time value into the  solution. Compare the answer to the amplitude calculated using the standard formula for converting the  solution into amplitude/phase angle form.   
c.  Assuming this is the model of a mass spring system, determine the speed of the mass as it passes  through equilibrium.   
d.  Convert the solution into the function g. Use g to plot the phase plane trajectory. What  type of curve is this trajectory?   
e.  Add to the trajectory the points corresponding to t = 0, 0.25, 0.5, 0.75, ... , 2.0.

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.    

a.   What is the pseudo-period of the oscillations?   
b.   What is the time constant?   
c.   Based upon your answer to part b estimate the time interval required for the oscillations to disappear  from view.   
d.   Plot the solution curve over the interval you named in part c.   
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?   
f.    Convert the solution into the function g. Use g to plot the phase plane trajectory.    
g.    Add to the trajectory the points corresponding to t = 0, 0.25, 0.5, 0.75, ... , 2.0.

Section 2. State Space

The forced oscillator is modeled with a non-autonomous equation. Solution trajectories are best viewed in state space.

1. Consider the driven IVP

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

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

a.    Plot the solution curve over the interval  0 ≤ t ≤ 120. What you witness in the plot is the phenomenon  called "beats". The output pulsates like this when the driver frequency is very close to the natural  frequency of the system.   
b.    Obtain the phase plane trajectory for this system. Use the same time interval.   
c.    Obtain the state space trajectory for this system. Use the same time interval.

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.     

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?   
b.    Plot the solution curve to verify your answer to part a.    
c.    Use h to obtain the phase plane and state space trajectories for this system.

3. Use PlotVectorField to obtain the direction field and a plot of the solution trajectory in Exercise 1 of Section 1 in Part IV. (The direction field can be drawn in phase space because the equation is  autonomous.)

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.

Section 3. Two Dimensional Systems

A two dimensional system defines a two dimensional vector field and a vector flow. The NDSolve  function outputs approximate solutions. Linear systems can be solved using standard  methods by reduction to one second order equation or by matrix methods (featured in the next section).

Load the PlotField package.

In[1]:=

<<Graphics`PlotField`

1. Obtain the phase portrait for the linear system

x' = -x  + 2 y
y' = 4 x + y    

Note. The phase portrait is a sketch of some solution curves and the direction field when the system is autonomous. Compare the picture you get to Figure 6.3.1 in Ledder (reproduced below).

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

Hint. Make the vector field using PlotVectorField. Use Table to generate lists of solutions satisfying initial conditions specified along the x and the y axes.

2. Continuing 1. Use DSolve to obtain the general solution to the system in Exercise 1. Call the solution  soln.

3. Add the two nullclines to the phase portrait in Exercise 1. (What is the stationary point?)  Hint. Go back and name the plot  PP.  Make the nullclines using ImplicitPlot (both can be  drawn at the same time). Call the nullcline plot NP and then use Show( PP, NP );

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.   
b.   Find the stationary point and add the nullclines to the direction field. Discuss the solutions based upon the picture you see.   
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

5. Use PlotVectorField to obtain the phase portrait for the system

x' = x - 2 y - 1
y' = x - y - 2

Add the stationary point and nullclines. Find the solution formulas using DSolve. Display the solutions in simplified form.
This is a linear system with periodic solutions (closed curves). Use the solution formulas to determine the period of the trajectories.

Section 4. Matrix Methods

Matrix and vector manipulation exercises provide practice via  eigenvalue and eigenvector calculations. The matrix exponential is a key tool for solving constant  coefficient linear systems of differential equations.

0. Pull down the Help menu, choose Help Browser... . Choose The Mathematica Book/Advanced Mathematics in Mathematica/Linear Algebra. Read the help pages.
Use Mathematica to do the following.

1. Enter the matrix (1   2   3)   4   5   6   7   8   9 with the name A.

a.    Find the determinant and the characteristic polynomial of A.   
b.   Enter the vector v = {2,3,5}. Calculate the vector w = Av using the entry w = A.v.
c.   Calculate the dot product of v and w using v.w. The product will be a scalar. Verify that the entry w.v yields the same scalar.  

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} } ]

a.   Find the characteristic polynomial of B and factor it.   
b.   Find B's eigenvalues with
Eigenvalues[B]. Name them lambda.
c.   Find B's eigenvectors using
Eigenvectors[B]. Name them V and define P to be the transpose of V.   
d.   Calculate
Inverse[P].B.P.

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.   
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)

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

v(t) = Mexp v(0)  

Section 5. The Laplace Transform

The Laplace transform and the inverse Laplace  transform are used to solve linear differential equations and systems. Piecewise defined functions are defined using the unit step function. Dirac delta functions in the driver can be handled via Laplace transforms.
Use Mathematica to do the following problems.

1. Obtain the Laplace transform of the following functions.

t cos(3t)    ,    e^(-t) sin(2t)    ,    ... p(t - π) (t + cos(t))    ,    sin(t) + DiracDelta(t - 3 π)

2. Obtain the inverse Laplace transform of the following functions.

s/(s^2 - 4)    ,     (s e^(-2s))/(s^2 - 4)    , & ...   ,    1/(s^2 + 4s - 4)    ,    e^(-3s)/s^4

3. Consider the following initial value problem    

y'' + y' + 9 y = UnitStep(t - 2 π )  ,  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.  
b.   Obtain the solution using the method of Laplace transforms. How does the solution formula compare to the  formula obtained in part a?

4. Consider the following initial value problem

y'' + y' + 9 y = UnitStep(t - 2 π ) - UnitStep(t - 5 π )  ,  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.  
b.   Obtain the solution using the method of Laplace transforms. How does the solution formula compare to the  formula obtained in part a?

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.  
b.   Obtain the solution using the method of Laplace transforms. How does the solution formula compare to the  formula obtained in part a?

6.* Use the the solution to 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.

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


Created by Mathematica  (December 6, 2004)