Section 5. The Laplace Transform

6.* Consider the following initial value problem

y'' + y' + 9 y = H(t - 2 Pi) - H(t - 5 Pi) + 2 Dirac(t - 9 Pi)  ,  y(0) = 0 , y'(0) = 1

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

> alias(H=Heaviside,delta=Dirac);

H, delta

> DE := diff(y(t),t,t) + diff(y(t),t) + 9*y(t) = H(t-2*Pi)-H(t-5*Pi)+2*delta(t-9*Pi);
inits := y(0)=0, D(y)(0)=1;

DE := (diff(y(t), `$`(t, 2)))+(diff(y(t), t))+9*y(t) = H(t-2*Pi)-H(t-5*Pi)+2*delta(t-9*Pi)

inits := y(0) = 0, D(y)(0) = 1

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

soln := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)-1/9*H(t-2*Pi)*exp(-1/2*t+Pi)*cos(1/2*35^(1/2)*t-35^(1/2)*Pi)-1/315*35^(1/2)*H(t-2*Pi)*exp(-1/2*t+Pi)*sin(1/2*35^(1/2)*t-35^(1/2)*Pi)+1/9*H(...soln := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)-1/9*H(t-2*Pi)*exp(-1/2*t+Pi)*cos(1/2*35^(1/2)*t-35^(1/2)*Pi)-1/315*35^(1/2)*H(t-2*Pi)*exp(-1/2*t+Pi)*sin(1/2*35^(1/2)*t-35^(1/2)*Pi)+1/9*H(...soln := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)-1/9*H(t-2*Pi)*exp(-1/2*t+Pi)*cos(1/2*35^(1/2)*t-35^(1/2)*Pi)-1/315*35^(1/2)*H(t-2*Pi)*exp(-1/2*t+Pi)*sin(1/2*35^(1/2)*t-35^(1/2)*Pi)+1/9*H(...soln := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)-1/9*H(t-2*Pi)*exp(-1/2*t+Pi)*cos(1/2*35^(1/2)*t-35^(1/2)*Pi)-1/315*35^(1/2)*H(t-2*Pi)*exp(-1/2*t+Pi)*sin(1/2*35^(1/2)*t-35^(1/2)*Pi)+1/9*H(...

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

[Plot]

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.

> plot( rhs(DE), t=0..40, title="Heaviside in Driver");

[Plot]

b. Obtain the solution using dsolve/method=laplace. How does the solution formula compare to the formula obtained in part a?

> soln2 := dsolve( {DE,inits}, y(t), method=laplace );

soln2 := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)+4/35*H(t-9*Pi)*35^(1/2)*exp(-1/2*t+9/2*Pi)*sin(1/2*35^(1/2)*(t-9*Pi))+1/630*(70-I*(-35*I+35^(1/2))*exp(-1/2*I*(-I+35^(1/2))*(t-2*Pi))+I*ex...soln2 := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)+4/35*H(t-9*Pi)*35^(1/2)*exp(-1/2*t+9/2*Pi)*sin(1/2*35^(1/2)*(t-9*Pi))+1/630*(70-I*(-35*I+35^(1/2))*exp(-1/2*I*(-I+35^(1/2))*(t-2*Pi))+I*ex...soln2 := y(t) = 2/35*exp(-1/2*t)*sin(1/2*35^(1/2)*t)*35^(1/2)+4/35*H(t-9*Pi)*35^(1/2)*exp(-1/2*t+9/2*Pi)*sin(1/2*35^(1/2)*(t-9*Pi))+1/630*(70-I*(-35*I+35^(1/2))*exp(-1/2*I*(-I+35^(1/2))*(t-2*Pi))+I*ex...

The Laplace solution is in complex form, probably because of a factorization of the characteristic polynomial over the complex field.

7.* Use unapply to convert the solution to 6 a into a function, g. Use g to plot the phase plane trajectory and as well as the state space trajectory for the system. Use the following options in the spacecurve procedure for the state space trajectory

axes=boxed, orientation=[-60,70], color=red, numpoints=600, labels=["Position","Velocity","Time"]

> g := unapply(rhs(soln),t):

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

[Plot]

> plots[spacecurve]( [g(t),D(g)(t),t,t=0..100],
                   title="State Space Trajectory",

                   axes=boxed, orientation=[-60,70], color=red,

                   numpoints=800,

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

[Plot]

8.* Continuing 7. Use the entry

plots[animate]( plot, [ [g(t),D(g)(t),t=0..T] ], T=0..40, frames=82);

to plot an animation of the phase space trajectory. Click on the plot and use the video controls in the context bar to play the animation and also to view it one frame at a time. Compare the trajectory movement to the time series for position.

Note. The slide control can also be used to control the display of the frames. Assuming t is measured in seconds, you can reduce the frame rate to 2 frames per second to view the trajectory in "real time".

> plots[animate]( plot, [ [g(t),D(g)(t),t=0..T] ], T=0..40, frames=82);

[Plot]