Part IV. Linear Differential Equations
Section 1. Linear Oscillators

Let's set the stage by recalling some of the things you are learning in your differential equations course.  A first order linear ordinary differential equation has the general form

a(t) y'(t) + b(t) y(t) = f(t)

or

a y' + b y = f

for short. Some equations like these were solved and analyzed in Part II. See the differential equation that was used to model Juanita's credit card payments and the equation modeling the temperature inside of a house on a  sunny day. A second order linear differential equation has the general form

a y'' + b y' + c y = f

The letters a, b, c, and f represent known functions of, say, time, and  y = y(t) is unknown. This is called a linear oscillator equation when a, b, and c are constant functions, a and c positive, b non-negative. Certainly the most famous example of a linear oscillator is the model for an object of mass m as it oscillates on a spring with spring constant k and possibly slowed down by a force proportional to its speed. The constant of proportionality is  called the damping constant and is often denoted by the letter b. Application of Newton's equally famous

Force = Mass x Acceleration

yields the equation

m y'' + b y' + k y = f

The function f is called the forcing function or the input because it usually represents an external force applied  to the mass. If f = 0, the equation is called "unforced" (or homogeneous in the general case). See Ledder, Chapter 3, Section 1.

Solve this too: DSolve and the unforced equation

Mathematica's DSolve function can be applied to an ordinary differential equation of any order. In particular it can  supply the solution formula for any unforced linear oscillator equation. Most of this section is devoted to  examples illustrating how to obtain solutions and plot their graphs.    

Example An object of mass 4.3 kilograms is attached to a spring with spring constant 15.2 Newtons per  meter. Find the general solution if the system is undamped.

In[162]:=

undampedDE = 4.3y''[t] + 15.2y[t] == 0
gensoln = DSolve[ undampedDE, y[t], t ]

Out[162]=

15.2 y[t] + 4.3 y^′′[t] 0

Out[163]=

{{y[t] ^(0. t) C[1] Cos[1.88013 t] + ^(0. t) C[2] Sin[1.88013 t]}}

Typical Problem Pull the mass 2 meters in the positive direction. Push it away from equilibrium with an initial speed of 3 meters per second. Plot the solution. How far does the mass move in the positive direction?  When does it pass through the equilibrium position for the first time? How fast is it going at that time?

This is an initial value problem, y(0) = 2, y'(0) = 3. We name the solution soln and use it to define the solution function, g. FullSimplify is applied to the solution formula to make it as simple as possible.

In[171]:=

soln = DSolve[ {undampedDE, y[0]==2, y'[0]==3}, y, t ];
g := soln[[1,1,2]]//FullSimplify
g[t]

Out[173]=

2. ^(0. t) Cos[1.88013 t] + 1.59564 ^(0. t) Sin[1.88013 t]

The graph of the solution is a shifted sine wave or, if you prefer, a shifted cosine wave. Such curves are referred to as sinusoids.

In[176]:=

Plot[ g[t], {t,0,6}, AspectRatio->1/3 ]

[Graphics:HTMLFiles/Math_P4S1_4.gif]

To find how far the mass moves in the positive direction we will find the first maximum value of  g for t positive. Name the t value tmax. The number ymax is obtained as g(tmax).

In[178]:=

tmax = t/.FindRoot[ g'[t]==0, {t,0.5} ]

Out[178]=

0.358172

In[180]:=

ymax = g[tmax]

General :: spell1 : Possible spelling error: new symbol name \"ymax\" is similar to existing symbol \"tmax\".  More…

Out[180]=

2.55853

So, the mass initially moves 0.558 meters in the positive direction to the position 2.558 meters from  equilibrium.  

To find when the mass passes through equilibrium, we solve the equation  g(t) = 0. The graph tells us that the first positive root of  g is near t = 1. We name it t0.

In[181]:=

t0 = t/.FindRoot[ g[t]==0, {t,1} ]

Out[181]=

1.19364

The velocity of the object at time t0 is g'(t0). It will be negative (see the graph), so apply the absolute value function Abs to get the speed.

In[182]:=

Abs[g'[t0]]*meterspersecond

Out[182]=

4.81036 meterspersecond

This information is plotted on the following graph. Position is red, velocity blue, and the dashed lines indicate the point where the mass passes through equilibrium for the first time.

In[190]:=

lines = ListPlot[ {{0,-4.81},{t0,-4.81},{t0,0}}, PlotJoined->True,
                  PlotStyle->Dashing[{0.02,0.02}] ];
Show[ lines,
      Plot[ {g[t],g'[t]}, {t,0,2},
            PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]} ] ,
      PlotRange->{-6,4}
    ]

[Graphics:HTMLFiles/Math_P4S1_10.gif]

The phase plane trajectory

The information contained in the last plot can also be displayed using what is known as the phase plane trajectory. This is a plot of the parametrized curve y = y(t) , v = y'(t) in two dimensional y,v space called the phase plane (v denotes velocity). Since g is defined so that y = g(t), the phase plane trajectory for the first 2  seconds is obtained by applying the ParametricPlot procedure to the list { g(t), g'(t) }. The setting AspectRatio->6/4 makes sure that the same scale is used on the two axes.

In[201]:=

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

[Graphics:HTMLFiles/Math_P4S1_11.gif]

The next plot shows the same trajectory, but with the phase points corresponding to position and velocity at  the 9 time values t = 0, 0.25, 0.5, ... , 2.0.

In[212]:=

Show[ traj, ListPlot[ Table[ {g[t],g'[t]}, {t,0,2,0.25}], PlotStyle->PointSize[0.03] ]]

[Graphics:HTMLFiles/Math_P4S1_12.gif]

The key to phase plane trajectories

Given a point on the trajectory, its first coordinate is the position  of the object and its second coordinate is its velocity. The nine points on the trajectory above correspond to the nine  pairs of points appearing on the position (red) and velocity (blue) curves displayed below. By a "pair" of points we mean two points corresponding to the same t value.

In[209]:=

pos = ListPlot[ Table[ {t,g[t]}, {t,0,2,0.25}],
                PlotStyle->PointSize[0.02] ];
vel = ListPlot[ Table[ {t,g'[t]}, {t,0,2,0.25}],
                PlotStyle->PointSize[0.02] ];
Show[ pos, vel,
      Plot[ {g[t],g'[t]}, {t,0,2},
            PlotStyle->{RGBColor[1,0,1],RGBColor[0,0,1]} ]
    ]

[Graphics:HTMLFiles/Math_P4S1_13.gif]

Simple Harmonic Motion

Simple harmonic motion is the name given to the motion of an undamped mass on a linear spring. The plot of  position versus time, called a "time series", is a sinusoid. The phase plane trajectory for simple harmonic  motion is an ellipse. Only a portion of the ellipse is shown in the trajectory plotted above. The time it takes a phase point to traverse the ellipse is the same time it takes the mass spring system to return to its original configuration. The motion is periodic. The extreme values on the time series for position correspond to the points on the phase plane  trajectory that are farthest from the vertical velocity axis. The extreme values on the time series for velocity  correspond to the points on the phase plane trajectory that are farthest from the horizontal position axis.

Underdamped, unforced motion

We now examine the geometry of unforced oscillations when the system is underdamped. Solutions still oscillate. The overdamped case is considered in the exercises. Consider the same mass and spring, but add a damping term with b = 1.0.

In[4]:=

dampedDE = 4.3y''[t] + y'[t] + 15.2y[t] == 0
gensoln = DSolve[ dampedDE, y[t], t]

Out[4]=

15.2 y[t] + y^′[t] + 4.3 y^′′[t] 0

Out[5]=

{{y[t] ^(-0.116279 t) C[2] Cos[1.87653 t] + ^(-0.116279 t) C[1] Sin[1.87653 t]}}

Observe that the general solution still contains a sine and a cosine, but now a decaying exponential term appears on each one. Factor it out and what remains is a sinusoid.

Typical Problem Pull the mass 2 meters in the positive direction and push it away from equilibrium with an  initial speed of 3 meters per second, just like before. Plot the solution and answer the same questions about its motion. Compare the answers to the ones obtained in the undamped case.

These are the same initial conditions. We call the solution function G.

In[6]:=

soln = DSolve[ {dampedDE, y[0]==2, y'[0]==3}, y, t];
G := soln[[1,1,2]]
G[t]

Out[8]=

^(-0.116279 t) (2. Cos[1.87653 t] + 1.72263 Sin[1.87653 t])

The graph:

In[225]:=

Plot[ G[t], {t,0,6}, AspectRatio->1/3]

[Graphics:HTMLFiles/Math_P4S1_17.gif]

This not a sinusoid because the oscillations are not of constant amplitude, they are decreasing. The time  constant for the exponential term is approximately 9 seconds (0.11 is about 1/9) so, over the 6 second time interval shown above, the oscillations decrease by a factor of about 6/9 x 1/3 = 2/9. The t value for the widest oscillation will  be called tmaxG.

In[227]:=

tmaxG = t/.FindRoot[ G'[t]==0, {t,0.5} ]
ymaxG = G[tmaxG]

Out[227]=

0.345926

Out[228]=

2.53067

The damped mass moves 0.531 meters in the positive direction (compared to 0.558 meters before) and it gets  there a tiny bit faster (the tmax value calculated earlier was 0.358). The time that the damped mass first passes through equilibrium will be denoted t0G.  

In[229]:=

t0G = t/.FindRoot[ G[t]==0, {t,1} ]

Out[229]=

1.21598

Recall that the undamped mass got to the equilibrium for the first time at  t = 1.194 seconds. It was traveling at  a speed of 4.809 meters a second. The next calculation shows that the damped mass not only took longer, but it is also going more slowly through equilibrium, traveling 4.300 meters per second.

In[230]:=

Abs[G'[t0G]]*meterspersecond

Out[230]=

4.30017 meterspersecond

This information goes on the the graph of the position and velocity, with dashed lines, just like before.

In[231]:=

lines = ListPlot[ {{0,-4.30},{t0G,-4.30},{t0,0}}, PlotJoined->True,
                  PlotStyle->Dashing[{0.02,0.02}] ];
Show[ lines,
      Plot[ {G[t],G'[t]}, {t,0,2},
            PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]} ] ,
      PlotRange->{-6,4}
    ]

[Graphics:HTMLFiles/Math_P4S1_23.gif]

The phase plane trajectory for underdamped oscillations

The motion of the damped system is no longer periodic and the phase plane trajectory is no longer an ellipse. It spirals inward. See the phase trajectory below where the time of motion has been extended to 6 seconds.

In[234]:=

traj = ParametricPlot[ {G[t],G'[t]}, {t,0,6}, PlotRange->{{-4,4},{-6,6}},
                AspectRatio->6/4 ]

[Graphics:HTMLFiles/Math_P4S1_24.gif]

The next picture adds the 25 points to the trajectory corresponding to position and velocity at time values   t =  0, 0.25, 0.5, ... , 6.0.

In[235]:=

Show[ traj, ListPlot[ Table[ {G[t],G'[t]}, {t,0,6,0.25}], PlotStyle->PointSize[0.03] ]]

[Graphics:HTMLFiles/Math_P4S1_25.gif]

And here are the corresponding pairs of points on the time series graphs.

In[12]:=

pos = ListPlot[ Table[ {t,G[t]}, {t,0,6,0.25}],
                PlotStyle->PointSize[0.01], DisplayFunction→Identity ];
vel = ListPlot[ Table[ {t,G'[t]}, {t,0,6,0.25}],
                PlotStyle->PointSize[0.01], DisplayFunction→Identity ];
Show[ pos, vel,
      Plot[ {G[t],G'[t]}, {t,0,6},
            PlotStyle->{RGBColor[1,0,1],RGBColor[0,0,1]},
            DisplayFunction→Identity ] ,
      AspectRatio→1/3, DisplayFunction→$DisplayFunction ]

[Graphics:HTMLFiles/Math_P4S1_26.gif]

Phase plane trajectories and vector fields

When the phase plane trajectory passes through the phase point (y, v) its velocity vector, which is tangent to the trajectory, is  < y', v' >. Since y' = v and  v' = y'', we have  

< y', v' > = < v, y'' > .

This tells us several things about the trajectories for a linear oscillator.   

1. All phase plane trajectories cross the y axis with a vertical tangent. This is because v = 0 on the y axis so the first component of the velocity vector is also 0.   

2. Similarly, trajectories move to the right when they are in the upper half plane, and move to the left in the  lower half plane.    

3. Since the differential equation determines y'' in terms of y and y' (= v), the trajectory's velocity vector can  be found without knowledge of the solution formula. Consider, for example, a second order linear constant coefficient homogeneous equation: (a not zero)

a y'' + b y' + c y = 0 .

Substitute y' = v and solve for y'' to obtain  

y'' = (-c y - b v)/a .  

Consequently, when a phase plane trajectory is at (y, v ), its velocity vector is

< y', v' > = < v , (-c y - b v)/a > .

This vector equation can be expressed equivalently as a system of two first order scalar differential equations.

y' = v                 
v' = (-c y - b v)/a

The system of differential equations is called autonomous because the two derivatives do not depend directly upon  the time variable. In this case, a plot of an array of tangent vectors in the y,v plane gives the same information  about the trajectories associated with a second order differential equation as a direction field gives for a first  order equation. Such an array is called the vector field associated with the differential equation.

Example. The vector field for the damped mass spring system in this section is defined below. We name it F.

In[247]:=

dampedDE

Out[247]=

15.2 y[t] + y^′[t] + 4.3 y^′′[t] 0

In[248]:=

F = { v, (-15.2y - v)/4.3 }

Out[248]=

{v, 0.232558 (-v - 15.2 y)}

Recall that Mathematica's PlotVectorField function will plot the vector field. It must be loaded into the kernel.

In[249]:=

<<Graphics`PlotField`

In[257]:=

vf = PlotVectorField[ F, {y,-5,5}, {v,-5,5}, Axes->True ]

[Graphics:HTMLFiles/Math_P4S1_29.gif]

The Show function will display this field and the last trajectory shown above.

In[261]:=

Show[ traj, vf, PlotRange->{{-4,4},{-6,6}} ]

[Graphics:HTMLFiles/Math_P4S1_30.gif]

Go to Next Section Previous Section


Created by Mathematica  (September 13, 2004)