Part IV. Linear Differential Equations
Section 5. The Laplace Transform
In this section, examples show how to use Maple's Laplace transform procedures to obtain transforms and their inverses. Applications to problems with discontinuous forcing functions are featured. See Ledder, Chapter 7.
The inttrans package
The standard integral transform procedures of applied mathematics are contained in Maple's inttrans package.
> | with(inttrans); |
We are interested in the procedures named laplace and invlaplace. The syntax to generate the Laplace transform of an expression f(t) as an expression F(s) is
laplace( f(t), t, s )
The inverse Laplace transform of F(s) to f(t) is accomplished with
invlaplace( F(s), s, t )
Four textbook examples:
1. The Laplace transform of
> | laplace( t*exp(2*t)*cos(3*t), t, s ); |
2. The inverse Laplace transform of
> | invlaplace( (s^2-1)/(s^3+8), s, t); |
3. Add this to your Laplace transform table: The transform of
> | laplace( t*exp(a*t)*sin(k*t), t, s); |
4. The transform of a discontinuous time function:
> | laplace( cos(t) + 3*t*2 - piecewise(t<3,0,1), t, s); |
"piecewise"? What is that all about?
Answer: The piecewise procedure makes functions that are piecewise continuous. The entry
piecewise(t<3,0,1)
defines the function that has the value 0 when t < 3 and the value 1 otherwise. It prettyprints using the "stacked" notation that is used in textbooks.
> | piecewise(t<3,0,1); |
See how it graphs.
> | plot( piecewise(t<3, 0, 1), t=0..6, -0.5..1.5); |
Piecewise versus Heaviside: Both are winners
If a piecewise function appears as part of an inverse Laplace transform, Maple uses the unit step function in the formula. This function is named Heaviside in honor of the man who was among the first to use it in engineering applications, Oliver Heaviside.
In Heaviside notation the function piecewise(t<3,0,1) is
Heaviside(t-3)
> | plot( Heaviside(t-3), t=0..6, -0.5..1.5); |
Maple's convert procedure can be used to convert form one form to the other.
> | convert(piecewise(t<3,0,1), Heaviside); |
> | convert(Heaviside(t-3), piecewise); |
Observe that Maple's version of the Heaviside function is undefined at t = 0. This can be changed if desired. For example, some applications require Heaviside(0) = 1. If this is the case, simply enter
> | Heaviside(0) := 1; |
Maple will make note of this value in a special table reserved for function values and then use it in all subsequent calculations. Compare the following output to the output we got to the same input above.
> | convert(Heaviside(t-3), piecewise); |
Of course, it makes no difference in the graphs (almost).
Warning: It does make a diffence when plotting points. Moreover, the fact that default Heaviside(0) is undefined can cause problems in the plotting of other graphs. Maple's plotter sometimes stops dead when it encounters "undefined". You will know it when you see it. Therefore, we will leave the value of Heaviside(0) at 1.
Be on the lookout for Heaviside when inverting an s-domain function containing expressions the form
, a positive.
> | invlaplace( s/(s^2+1)+6/s^2-exp(-3*s)/s, s, t); |
And, be prepared to use Heaviside to deal with piecewise continuous functions. It happens all the time in engineering and science applications.
Typical driver: Off, then on, then off again
Suppose, for example, an undamped mass spring system is set into motion from equilibrium with velocity 3 meters per second, then driven with a sinusoid for the 10 second interval between t = 10 and t = 20. See the graph below.
> | f := t -> (Heaviside(t-10) - Heaviside(t-20))*20*sin(Pi/4*t);
plot( f(t), t=0..40, -22..22); |
Here is the Laplace transform of the driver.
> | laplace(f(t), t, s); |
For the forced differential equation we will use m = 40 (kilograms), b = 0, k = 10.
> | DE := 40*diff(y(t),t,t) + 10*y(t) = f(t); |
The dsolve solution to the IVP is obtained first.
> | soln := dsolve( {DE, y(0)=0, D(y)(0)=3}, y(t) ); |
dsolve did not use the method of Laplace transforms to solve this equation.
Question: What method did dsolve use?
Answer: dsolve sees that the equation has constant coefficients, solves the homogeneous equation, then uses an integral formula derived from the method of variation of parameters. See Ledder, Chapter 4, Section 6, especially Example 2 and Exercise 16.
Question: How do you know?
Answer: I set the "infolevel" for dsolve to level 3, like this. (The default infolevel for any procedure is 1.)
> | infolevel[dsolve] := 3; |
Then I ran the dsolve procedure again, suppressing output, and got the following information.
> | soln := dsolve( {DE, y(0)=0, D(y)(0)=3}, y(t) ): |
`Methods for second order ODEs:`
`--- Trying classification methods ---`
`trying a quadrature`
`trying high order exact linear fully integrable`
`trying differential order: 2; linear nonhomogeneous with symmetry [0,1]\
`
`trying a double symmetry of the form [xi=0, eta=F(x)]`
`Try solving first the homogeneous part of the ODE`
` -> Tackling the linear ODE "as given":`
` checking if the LODE has constant coefficients`
` <- constant coefficients successful`
` <- successful solving of the linear ODE "as given"`
` -> Determining now a particular solution to the non-homogeneous ODE`\
` building a particular solution using variation of parameters`
` <- solving first the homogeneous part of the ODE successful`
That's how I know what dsolve did. And that is how you can discover some of the details about most of Maple's procedures.
Before continuing we set the infolevel for dsolve back to level 1. Too much information is not always useful.
> | infolevel[dsolve] := 1: |
Question: Well then, can Maple apply the method of Laplace transforms using dsolve?
Answer: Sure, simply insert the unknown function right after {DE,inits}, and add the equation
method=laplace
> | laplace_soln := dsolve( {DE, y(0)=0, D(y)(0)=3}, y(t), method=laplace ); |
To see the effect of using the Laplace transform method, compare this form of the solution to the form obtained earlier. You might like to convert the first solution formula to piecewise form for easy interpretation.
> | convert(soln,piecewise); |
The graph of the solution shows that the input between 10 and 20 (shown scaled down to fit and in blue) caused the amplitude of the oscillations to decrease from 6 to 4 meters (approximately). The phase of the oscillations has also been affected, but not the period.
> | plot([rhs(soln),f(t)/10], t=0..40, color=[red,blue],
title="Blue curve is the scaled input"); |
Use an alias
After a while typing the word Heaviside all the time gets to be a chore. Why can't we just use an H? Well, we can. The following entry tells Maple that from now on (in this session) every time we type H what we mean is Heaviside.
> | alias(H=Heaviside); |
Maple will also use H for Heaviside in the display of its outputs.
> | laplace_soln; |
The Dirac delta: Dirac
The Dirac delta function is written as Dirac in Maple.
> | laplace(Dirac(t-3),t,s); |
> | invlaplace( 1, t, s); |
You might like to use a delta, like your textbook probably does. No problem, make an alias.
> | alias(delta=Dirac); |
The output is a sequence listing all of the aliases in effect at the present time in the session.
Lets drive the mass spring system introduced above with a blow of magnitude 200 in the positive direction at time t = 12.
> | DE2 := lhs(DE) = 200*delta(t-12); |
The solution formula, using the same initial conditions and the method of Laplace transforms, looks like this.
> | lap_soln2 := dsolve( {DE2, y(0)=0,D(y)(0)=3}, y(t), method=laplace); |
The following plot displays the solution curve (output) and a scaled arrow placed at the time of the application of the blow (input). See the plots[arrow] help page to learn about arrows in plots.
> | with(plots):
input := arrow([12,0],[0,6],width=0.4, head_width=1, color=blue): output := plot(rhs(lap_soln2), t=0..40): display( input, output, title="A blow at t = 12"); |
The effect of the blow is to cause a kink in the curve for position as a function of time. There will also be a jump discontinuity on the curve for velocity.
The phase plane plot shows the discontinuity in velocity as a vertical line joining the initial elliptical trajectory to the elliptical trajectory describing the motion of the system after the blow.
> | plot( [rhs(lap_soln2),diff(rhs(lap_soln2),t),t=0..22],
view=[-20..20,-10..10], axes=framed, labels=["POSITION","VELOCITY"], labeldirections=[horizontal,vertical]); |
Applying the Laplace transform method to a linear system
The equation method=laplace can also be added to dsolve when it is applied to a system of linear equations (constant coefficients), driven or undriven. However, if the system is driven by discontinuous forcing functions, the Laplace transform solution formulas may contain complex exponentials.
Consider the following IVP
x' = x - y + 10 H(t-2)
y' = x + y
x(0) = 1, y(0) = 1
> | sys := diff(x(t),t) = x(t) - y(t) + 10*H(t-2),
diff(y(t),t) = x(t) + y(t); inits := x(0)=1, y(0)=1; |
Maple's Laplace transform solution outputs in complex form.
> | lap_soln := dsolve( {sys,inits}, [x(t),y(t)], method=laplace); |
The solution obtained using the default variation of parameters technique does not.
> | soln := dsolve( {sys,inits}, [x(t),y(t)]); |
Plot the solution trajectory
The following entry shows how to get the "soln" output into trajectory form for plotting.
> | v := subs(soln,<x(t),y(t)>); |
The trajectory v is plotted below along with a point plotted for every quarter unit of time. See if you can identify the point where the Heaviside driver starts to have an effect on the output.
> | Pts := [eval([v[1],v[2]],t=k/4)$k=0..16]:
plot( [[v[1],v[2],t=0..4],Pts] , x=-15..25, y=-25..5, style=[line,point], color=[green,black], axes=framed); |
Go to: Previous Section