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.

> DE := diff(y(t),t,t) + diff(y(t),t) + 4*y(t) = 0;
inits := y(0)=2, D(y)(0)=-3;

DE := (diff(y(t), `$`(t, 2)))+(diff(y(t), t))+4*y(t) = 0

inits := y(0) = 2, D(y)(0) = -3

> soln := dsolve( {DE,inits} );

soln := y(t) = -4/15*15^(1/2)*exp(-1/2*t)*sin(1/2*15^(1/2)*t)+2*exp(-1/2*t)*cos(1/2*15^(1/2)*t)

The system is underdamped.

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

The pseudo-period is

> 2*Pi/(sqrt(15)/2): evalf[4](%)*time_units;

3.246*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.

The oscillations will disappear in about 5 x 2 = 10 time units.

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

> plot( rhs(soln), t=0..10);

[Plot]

e. Add to the curve in part d the curves defined by Aexp(-t/2) and -Aexp(t/2) where

A = sqrt(4+16/15) .

Make them blue. What is the significance of these curves? Where did the formula for A come from?

> A := sqrt(4 + 16/15):
plot( [rhs(soln),A*exp(-t/2),-A*exp(-t/2)], t=0..10, color=[red,blue$2]);

[Plot]

The blue curves are the envelopes for the oscillations.

The value of A is derived from the standard formula used to convert the sum of a sine and a cosine into a sinusoid.

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

> g := unapply(rhs(soln),t):
plot( [g(t),D(g)(t),t=0..10] );

[Plot]

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

> Points := [g(t/4),D(g)(t/4)] $ t=0..8:
plot( [ [g(t),D(g)(t),t=0..10], [ Points ] ], style=[line,point],

       color=[red,black] );

[Plot]

4.* Use DEplot to draw the direction field in phase space for the damped system. Then add the solution trajectory.

> DEsystem := diff(y(t),t) = v(t), diff(v(t),t) = -v(t) - 4*y(t):

> DEtools[DEplot]( {DEsystem}, [y(t),v(t)], t=0..10, y=-2..2, v=-4..2,
                 dirgrid=[11,11], arrows=slim);

[Plot]

> DEtools[DEplot]( {DEsystem}, [y(t),v(t)], t=0..10, y=-2..2, v=-4..2,
                 dirgrid=[11,11], arrows=slim, {[y(0)=2,v(0)=-3]},

                 linecolor=blue, stepsize=0.1);

[Plot]

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 g using unapply. Suppress the output for soln and the definition of g and then enter

evalf[3]( 'g(t)' = g(t) )

to see a nice looking representation of the solution formula.

> DE := diff(y(t),t,t) + 0.1*diff(y(t),t) + 4*y(t) = cos(1.8*t);
inits := y(0)=2, D(y)(0)=-3;

DE := (diff(y(t), `$`(t, 2)))+.1*(diff(y(t), t))+4*y(t) = cos(1.8*t)

inits := y(0) = 2, D(y)(0) = -3

> soln := dsolve( {DE,inits} ):
g := unapply(rhs(soln),t):

evalf[3]('g(t)'=g(t));

g(t) = -1.75*exp(-0.500e-1*t)*sin(2.00*t)+.754*exp(-0.500e-1*t)*cos(2.00*t)+1.25*cos(1.80*t)+.295*sin(1.80*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 is 1/0.05 time units.

> 1/0.05*time_units;

20.00000000*time_units

It will take approximately 100 seconds for the beats to disappear from the solution curve.

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

> plot( g(t), t=0..100 );

[Plot]

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

> plot( [g(t),D(g)(t),t=0..100], title="Phase Plane");

[Plot]

> plots[spacecurve]( [g(t),D(g)(t),t,t=0..100], numpoints=600,
     orientation=[-60,70], axes=boxed,

     labels=["Position","Velocity","Time"],

     color=red, title="State Space" );

[Plot]

5.* Apply DEsystem it to the differential equation in Exercise 2 in this section and use DEplot to obtain the solution trajectory.

> DEsystem := diff(y(t),t) = v(t),
           diff(v(t),t) = -0.1*v(t) - 4*y(t) + cos(1.8*t):

DEtools[DEplot]( {DEsystem}, [y(t),v(t)], t=0..100, y=-3..3, v=-6..5,

                 linecolor=blue, stepsize=0.1, {[y(0)=3,v(0)=-2]} ,

                 title="Phase Space");

[Plot]

6.* Continuing 5. Apply DEplot3d to obtain the state space trajectory for the IVP in Exercise 2.

> DEtools[DEplot3d]( {DEsystem}, [y(t),v(t)], t=0..100, y=-3..3,
                 v=-6..5, scene=[y,v,t], orientation=[-60,70],

                 axes=boxed, labels=["Position","Velocity","Time"],

                 title="State Space" ,

                 linecolor=blue, stepsize=0.1, {[y(0)=3,v(0)=-2]} );

[Plot]

Section 3. Two Dimensional Systems

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 DEplot to do the following.

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

> with(DEtools):

> DEsystem := diff(x(t),t) = x(t)*(1 - y(t) - 1.9*x(t)),
           diff(y(t),t) = y(t)*(1 - x(t) - 1.5*y(t));

DEsystem := diff(x(t), t) = x(t)*(1-y(t)-1.9*x(t)), diff(y(t), t) = y(t)*(1-x(t)-1.5*y(t))

> DEplot( {DEsystem}, [x(t),y(t)], t=0..10, x=0..1, y=0..1,
                 arrows=slim, dirgrid=[13,13]);

DF := %:

[Plot]

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

> stpt := solve( {1-y-1.9*x,1-x-1.5*y},{x,y}):
StatPoint := subs(stpt,[x,y]);

Nullclines := plots[implicitplot]( [1-y-1.9*x,1-x-1.5*y], x=0..1,

                              y=0..1, color=green, thickness=2):

StatPoint := [.2702702703, .4864864865]

> plots[display]( DF, Nullclines, plot( [StatPoint], style=point,
                                      color=black));

[Plot]

Solutions will stream towards the stationary point, moving as indicated in the four regions determined by the nullclines.

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.2*cos(Pi/6*k),y(0)=0.5+0.2*sin(Pi*k/6)] $ k=0..11 }

> DF := DEplot( {DEsystem}, [x(t),y(t)], t=-10..10, x=0..1, y=0..1,
   arrows=slim, dirgrid=[13,13], linecolor=blue,

   { [x(0)=0.3+0.2*cos(Pi/6*k),y(0)=0.5+0.2*sin(Pi*k/6)]$k=0..11 }):

plots[display]( DF, Nullclines, plot( [StatPoint], style=point,

               color=black));

[Plot]

Section 4. Matrix Methods

Section 5. The Laplace Transform