Part IV. Linear Differential Equations
Section 5. The Laplace Transform
In this section, examples show how to use Mathematica to obtain Laplace transforms and their inverses. The transform method for solving an ode and a linear system is illustrated. Applications to problems with discontinuous forcing functions are featured. See Ledder, Chapter 7.
Laplace Transforms and the Inverse Laplace Transform
We are interested in the Mathematica functions named LaplaceTransform and InverseLaplaceTransform. The syntax to generate the Laplace transform of an expression f(t) as an expression F(s) is
LaplaceTransform[ f(t), t, s ]
The inverse Laplace transform of F(s) to f(t) is accomplished with
InverseLaplaceTransform[ F(s), s, t ]
Four textbook examples:
1. The Laplace transform of t
In[175]:=
LaplaceTransform[t*Exp[2*t]*Cos[3*t],t,s]
Out[175]=
2. The inverse Laplace transform of .
In[177]:=
InverseLaplaceTransform[ (s^2-1)/(s^3+8),s,t]
Out[177]=
3. Add this to your Laplace transform table: The transform of t .
In[181]:=
LaplaceTransform[ t*Exp[a*t]*Sin[k*t],t,s]
Out[181]=
4. The transform of a discontinuous time function: cos(t) -
In[194]:=
LaplaceTransform[ Cos[t] - t - UnitStep[t-3],t,s]
Out[194]=
"UnitStep"? What is that all about.
Answer: The UnitStep function has its values defined as follows:
UnitStep[t] = 0 when t is negative
UnitStep[t] = 1 when t is positive.
See its graph below.
In[191]:=
Plot[ UnitStep[t], {t,-2,2}, PlotRange->{-0.5,1.5},
PlotStyle->Thickness[0.015], AspectRatio->2/4 ]
Consequently, UnitStep[t-3] is zero for t < 3 and 1 for t > 3. The UnitStep function and functions defined using it (like the one whose Laplace transform was found in 4) are called "piecewise continuous". The time function in Example 4 has a jump discontinuity at t = 3. Its graph looks like this.
In[199]:=
Plot[ Cos[t] + t - UnitStep[t-3], {t,-1,5}, PlotRange->{-1,3} ]
Be on the lookout for UnitStep when inverting an s-domain function containing expressions of the form
In[202]:=
InverseLaplaceTransform[ s/(s^2+1) - 1/s^2 - Exp[-3*s]/s, s, t]
Out[202]=
And, be prepared to use UnitStep 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, that 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.
In[45]:=
f[t_] := 20*(UnitStep[t-10] - UnitStep[t-20])*Sin[Pi*t/4]
Plot[ f[t], {t,0,40}, PlotRange->{-22,22}, AspectRatio->1/6,
PlotStyle->Thickness[0.005]]
Here is the Laplace transform of the driver.
In[3]:=
LaplaceTransform[f[t],t,s]
Out[3]=
For the forced differential equation we will use m = 40 (kilograms), b = 0, k = 10.
In[1]:=
DE = 40*y''[t] + 10*y[t] == f[t]
Out[1]=
The DSolve solution is messy.
In[20]:=
soln = DSolve[ {DE, y[0]==0, y'[0]==3}, y[t], t ]
Out[20]=
But it simplifies considerably.
In[22]:=
soln = Simplify[soln]
Out[22]=
The Laplace transform solution is much nicer looking. We can obtain it as follows.
Transform DE to the s-domain
In[32]:=
LaplaceTransform[DE,t,s]
Out[32]=
Substitute the initial values, and let Y replace LaplaceTransform[y[t],t,s]
In[33]:=
%/.{y[0]->0,y'[0]->3,LaplaceTransform[y[t],t,s]->Y}
Out[33]=
Solve for Y
In[34]:=
Solve[%,Y]
Out[34]=
Invert to get the solution formula (Compare to the simplified solution formula displayed above.)
In[35]:=
InverseLaplaceTransform[Y/.%,s,t]
Out[35]=
Here is a plot of the solution and the driver (between 10 and 20). The input graph has been scaled down to fit nicely into the plot. The driver causes 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.
In[44]:=
Plot[ {y[t]/.soln,f[t]/7}, {t,0,40}, AspectRatio→1/3]
The Dirac delta: DiracDelta
The Dirac delta function is written as DiracDelta in Mathematica.
In[47]:=
LaplaceTransform[DiracDelta[t-3], t, s]
Out[47]=
In[48]:=
InverseLaplaceTransform[1, s, t]
Out[48]=
Lets drive the mass spring system introduced above with a blow of magnitude 200 in the positive direction at time t = 12. Note that DE[[1]] can be used to refer to the left side the DE equation.
In[2]:=
DE2 = DE[[1]] == 200*DiracDelta[t-12]
Out[2]=
The solution formula, using the same initial conditions, looks like this.
In[3]:=
soln = DSolve[ {DE2, y[0]==0, y'[0]==3}, y[t], t]
Out[3]=
The following plot displays the solution curve (output) and vertical line placed at the time of the application of the blow (input).
In[6]:=
Show[ Plot[ y[t]/.soln, {t,0,40} ],
ListPlot[ {{12,0},{12,5}}, PlotJoined->True,
PlotStyle->Thickness[0.015] ], AspectRatio→1/3 ]
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 below 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.
In[113]:=
Trajectory = {y[t]/.soln[[1]],D[y[t]/.soln[[1]],t]};
ParametricPlot[ Trajectory, {t,0,22}, AxesOrigin->{-20,-10},
PlotRange->{{-20,20},{-10,10}},
AxesLabel->{"POSITION","VELOCITY"}]
Applying the Laplace transform method to a linear system
The Laplace transform method can also be 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 UnitStep(t - 2)
y' = x + y
x(0) = 1, y(0) = 1
In[23]:=
DEsystem = {x'[t] == x[t] - y[t] + 10*UnitStep[t-2],
y'[t] == x[t] + y[t] }
Out[23]=
Here is the DSolve solution for the initial conditions x(0) = 1 and y(0) = 1.
In[24]:=
DSsoln = DSolve[ Join[DEsystem,{x[0]==1, y[0]==1}], {x[t],y[t]}, t]
Out[24]=
Is simplifies as follows.
In[25]:=
DSsoln = Simplify[{x[t],y[t]}/.DSsoln]
Out[25]=
The Laplace transform solution will arrive in complex form. See below.
Transform the system
In[40]:=
LaplaceTransform[DEsystem, t, s ]
Out[40]=
Substitute the initial values and X and Y for the transforms of the solution
In[41]:=
%/.{x[0]->1,y[0]->1,LaplaceTransform[x[t],t,s]->X,LaplaceTransform[y[t],t,s]->Y}
Out[41]=
Solve for X and Y
In[42]:=
Tsoln = Solve[%,{X,Y}]
Out[42]=
Invert to obtain the solution trajectory (Note that the solution is in terms of complex exponentials)
In[43]:=
Lsoln = InverseLaplaceTransform[{X,Y}/.%, s, t]
Out[43]=
The function called ComplexExpand will change the complex exponentials to Cartesian form, but the output is very long. We suppress it, then Simplify. Compare the output to the simplified DSolve solution shown above.
In[44]:=
ComplexExpand[Lsoln];
Simplify[%]
Out[45]=
Plot of the solution trajectory
The solution trajectory is plotted below along with a point plotted for every quarter unit of time. See if you can identify the point where the UnitStep driver starts to have an effect on the output.
In[68]:=
Show[ ListPlot[Table[DSsoln[[1]],{t,0,4,1/4}],PlotStyle->PointSize[0.02]],
ParametricPlot[ DSsoln, {t,0,4}, PlotStyle->RGBColor[0,1,0]],
AspectRatio->1/1, PlotRange->{{-15,25},{-25,5}},
AxesOrigin->{-15,-25}, AxesLabel->{"x","y"} ]
Go to Previous Section
Created by Mathematica (November 25, 2004)