7-5-S-Laplace Transform Issues and Opportunities in Mathematica

The Laplace Transform is a mathematical construct that has proven very useful in both solving and understanding differential equations. We introduce it and show its power here.

Definition of Laplace Transform of a function f(t) if the integral converges:

ℒ(f(t)) (s)= F(s) = (∫_0^∞)(e^(-s t))f(t) dt

In order to evaluate such an integral we might be using integration by parts for f(t) of a different "sort" than exponential functions.  However, we turn to Mathematica to do the integrations for us and use its own powerful computation, look-up, and pattern matching algorithms to do Laplace Transforms.

This "strange" integral transforms most functions (if the integral converges) from a t domain to an s domain where s is (in our case) a real positive number.

Integrate[Exp[- s t]f[t],{t,0, Infinity}]

Integrate[Exp[- s t]Cos[3t],{t,0, Infinity}]

Notice the restriction that Re[s]>0 in order for us to have convergence to the quantity, s/(9+s^2). Thus s/(9+s^2) IS the Laplace Transform of the function cos(3t).

Tables of these exist and we could visit them, however since Mathematica can do these instantly without us looking them up we  do our computations of Laplace Transforms in Mathematica.

Often the s domain is called the frequency domain, and the t domain is called the time domain, for if the following integral - the Laplace Transform - is to make sense with respect to units then the units of s must be 1/Time, for in this way s*t has no units and so Exp[-s*t] has no units. Thus the Laplace Transform integral is obtained by just integrating a function of t over an interval of time,  [0, ∞) - resulting in a function of s, we call it ℒ (f(t))(s) = F(s).

F[s_] = Integrate[Exp[- s t]f[t],{t,0, Infinity}]

Mathematica can compute these Laplace Transforms and we see how to do this.  Here are the definitions brought up by query from Mathematica.  If one needs more, such as an example, then go to Help Menu and type "LaplaceTransform" and you will have a window offering up a definition with  an example worked out at the bottom.OR just type one of the following

??LaplaceTransform

We have a tougher trick in undoing the Laplace Transform and this is done using the InverseLaplaceTransform command.

??InverseLaplaceTransform

So here we go . . . taking Laplace Transforms of functions . . . . We Build Up a Table of Laplace Transforms and attempt to ascertain patterns.

LaplaceTransform[1,t,s]

LaplaceTransform[t,t,s]

LaplaceTransform[t^2,t,s]

LaplaceTransform[t^3,t,s]

LaplaceTransform[t^4,t,s]

What do you believe we might get for the Laplace Transform of (t^n)? Of course, we could stop and offer a formal induction proof of  general formula, but we proceed to take the Laplace Transform of other classes of functions.

Trig functions

LaplaceTransform[Sin[3t],t,s]

LaplaceTransform[Sin[5t],t,s]

LaplaceTransform[Cos[3t],t,s]

LaplaceTransform[Cos[5t],t,s]

See a pattern here?

. . . and exponential functions.

LaplaceTransform[Exp[3 t],t,s]

LaplaceTransform[Exp[5 t],t,s]

LaplaceTransform[Exp[-7 t],t,s]

Incidentally, the transform of functions like f(t) = (e^(a t)) and f(t) = k are easy by hand efforts! Try one.

Continuing in this manner of discovery we could fill in a typical Laplace Transform Table. A sample, but short and simple form, table is the following:

We could combine trig and exponential functions . .

LaplaceTransform[Exp[5 t] Sin[3 t],t,s]

. . . and in general to match the table entry.

LaplaceTransform[Exp[a t] Sin[b t],t,s]

If we are given the Laplace Transform Y(s) how might we obtain the original, i.e. the y(t) whose Laplace Transform is Y(s). We try it on some of our observed Transforms and some others functions.

InverseLaplaceTransform[1/(s- 4),s,t]

InverseLaplaceTransform[s-3/((s-3)^2 + 4),s,t]

InverseLaplaceTransform[1/s^7,s,t]

InverseLaplaceTransform[s/(s^2 + 4)+ 1/(s- 4),s,t]

InverseLaplaceTransform[s/(s^2 - 4),s,t]

ysol[t_] = InverseLaplaceTransform[(s+3)/((s+3)^2 + 4),s,t]

Plot[ysol[t],{t,0,2}]

InverseLaplaceTransform[(s+3)/((s+3)^2 + 4),s,t]

Here we need to use Euler' s formula, (e^(i θ)) = cos(θ) + i sin(θ), invoked by Mathematica’s Complex Expand command to offer a “real” presentation of our transform.

ComplexExpand[InverseLaplaceTransform[(s+3)/((s+3)^2 + 4),s,t]]

Further, we invoke Mathematica's TrigReduce command to use trigonometric identities to reduce our Inverse Laplace Transform to a simple form from which we could see so much, in this case a damped oscillator!

TrigReduce[ComplexExpand[InverseLaplaceTransform[(s+3)/((s+3)^2 + 4),s,t]]]

Let us take the Laplace Transform of general functions  and derivatives . . . .

LaplaceTransform[y[t],t,s]

LaplaceTransform[y'[t],t,s]

We take Laplace Transform of second derivative . . .

LaplaceTransform[y''[t],t,s]

We take Laplace Transform of third derivative . . .

LaplaceTransform[y'''[t],t,s]

We take Laplace Transform of fourth derivative . . .

LaplaceTransform[y''''[t],t,s]

Let us give LaplaceTransform[y[t],t,s] a name, say Y[s]. Then things look nicer.

LaplaceTransform[y'[t],t,s]/.{LaplaceTransform[y[t],t,s]->Y[s]}

LaplaceTransform[y''[t],t,s]/.{LaplaceTransform[y[t],t,s]->Y[s]}

Can we take the Laplace Transform of a whole  differential equations? Let's try.

eq = LaplaceTransform[a y''[t] + b y'[t]  + c y[t]==Sin[t],t,s]

Let us try to replace LaplaceTransform[y[t],t,s] with Y[s] here too.

eq= LaplaceTransform[a y''[t] + b y'[t]  + c y[t]==Sin[t],t,s]/.{LaplaceTransform[y[t],t,s]->Y[s]}

We can solve for the Laplace of the function in the above differential equation.

YSol[s_] =Y[s]/.Solve[eq,Y[s]][[1]]

Now if we could compute the Inverse Laplace Transform of Y (s) we would have the solution to our differential equation, y (t) . Let us do this. Our solution appears a bit on the “wild side” involving the constants a, b, and c, and the all too familiar radical(Sqrt[(b^2) - 4ac]) one obtains from the characteristic equation of the homogeneous portion of our original differential equation a y’’(t) + b y’(t)   + c y(t)  = sin(t)  in using eigenvalue solution strategies.

InverseLaplaceTransform[YSol[s],s,t]//Simplify

The Transfer function of a differential equation such as a*y’’(t) + b*y’(t)   + c*y(t)  = sin(t)  is Q(s) = 1/ (a*(s^2) + b*s + c)  (1 over the denominator in the Laplace Transform of the solution)  and notice the coefficients a, b, and c, of the LHS or of the system - not the driver sin(t). This characterizes the system response and the roots of the characteristic equation, a*(s^2) + b*s + c = 0, are the characteristic roots or eigenvalues of this system. Engineers get to "know" systems and the systems' behaviors by their Transfer functions and the eigenvalues; electrical engineers  call the eigenvalues "poles."

Indeed, if we let y(0) = 0 and y'(0) = 0 then YSol[s] = L[s]*Q[s], and so Q[s] = Ysol[s]/L[s], where YSol[s] is the output's Laplace Transform and L[s] is the input's Laplace Transform. This effectively shows that the Laplace Transform of the Solution is “simply” a multiplication of Laplace Transform of the input or driver function, i.e. L[s], and the transfer function Q[s] = Q(s) = 1/ (a*(s^2) + b*s + c) . The latter contains all the “action” of the differential equation model terms   a*y’’(t) + b*y’(t)   + c y(t) while the former, L[s], contains the “action” of the driver function, in this case sin(t).

So when we actually use reasonable initial conditions (y(0) = 0, no initial displacement, and y’(0) = 0, no initial velocity - depending upon the driver to get some activity going) on the inverse Laplace Transform, YSol[s], of our differential equation a*y’’(t) + b*y’(t)   + c*y(t)  = sin(t)

YSol[s]/.{y[0]->0,y'[0]->0}

We try a specific example, one we know we can solve in several other ways - by hand, with Mathematica DSolve.  We consider 

		y''(t) + 6 y'(t) + 5 y(t) = sin(t),    y(0) = 4,    y'(0) = 0.

First we use DSolve and some formatting to grab the solution for comparison.

ysol[t_] = y[t]/.DSolve[{1 y''[t] + 6 y'[t]  + 5 y[t]==Sin[t],y[0]==4,y'[0]==0},y[t],t][[1]]//Expand

We consider  this same differential equation and apply Laplace Transforms to all terms, i.e. both sides of the differential equation.

eq = LaplaceTransform[1 y''[t] + 6 y'[t]  + 5 y[t]==Sin[t],t,s]/.{LaplaceTransform[y[t],t,s]->Y[s],y[0]->4,y'[0]->0}

Now we solve for Y[s], the Laplace Transform of y[t] in order to get ready for the inverse Laplace Transform.

L[s_] =Y[s]/.Solve[eq,Y[s]][[1]]

Now we apply the inverse transform to get the solution for y[t].

l[t_] = InverseLaplaceTransform[L[s],s,t]//Expand

Notice the two portions of this soloution, namely the Transient solution -(105/104) (e^(-5 t))+((41 (e^(-t)))/8)and the Teady State -((3 Cos[t])/26)+((Sin[t])/13).  These would be produced by various hand techniques or Mathematica’s DSolve approach which we show here.

How does that compare to the DSolve solution? Exactly the same as we see when, again, we use DSolve to obtain a solution.

ysol[t_] = y[t]/.DSolve[{1 y''[t] + 6 y'[t]  + 5 y[t]==Sin[t],y[0]==4,y'[0]==0},y[t],t][[1]]//Expand

So what is the big deal?  Who would someone choose Laplace Transforms to solve a problem if DSolve or by hand can do it directly? 

Here is the point (well one of the many points of advantage) of Laplace Transforms. They convert calculus to algebra, then we play in algebra land, and then we return to calculus land. More formally, the Laplace Transform transforms the problem from the time domain (t) to the frequency domain in (s).

We consider two additional "Real Good" reasons for Laplace Transforms

Unit Step or Heaviside Function - what does this function do.

f1[x_] = UnitStep[x]

Plot[f1[x],{x,-8,8},PlotStyle->{Thickness[.01]}]

f2[x_] = UnitStep[x-3]

Plot[f2[x],{x,-8,8},PlotStyle->{Thickness[.01]}]

We could define our own unit step function with conditional or if statements as follows, but Mathematica already has it defined so we accept UnitStep[x] in Mathematica.
      
      unitstep(x) = 1, if x>0, and 0 otherwise.

unitstep[x_] = If[x>0,1,0]

Plot[unitstep[x],{x,-8,8},PlotStyle->{Thickness[.01]}]

Dirac Delta or Impulse Function

g[x_] = DiracDelta[x]

Strange plot . . . .

Plot[g[x],{x,-4,4},PlotStyle->{Thickness[.01]}]

We examine  the Dirac Delta function some more.  This function represents an instantaneous change.  E.g., we change the salt concentration by "dumping" in more salt ALL AT ONCE!  OR we "shock" a spring mass giving it an instantaneous force ALL AT ONCE!

The integral of an interval containing x = 0 for the DiracDelta[x] function is 1. What this says is the impulse offered by a driver of the form DiracDelta[x] imparts force of size 1 unit to our system, BUT all at once like the quick bang of a hammer, rather than a gradual application of force.

We check out some properties of the DiracDelta[x] function.

Integrate[g[x],{x,-1,1}]

Integrate[g[x],{x,.5,1}]

g[.1]

g[0]

g[2]

Transforms of UnitStep and Dirac Function - where does this fit into our Table of Laplace transforms?

f[t_] = UnitStep[t];

LaplaceTransform[f[t],t,s]

g[t_] = DiracDelta[t];

LaplaceTransform[g[t],t,s]

UnitStep function represents a "step" or finite immediate change in something. DiracDelta function represents an instantaneous action at one point. . . zero elsewhere.

del[x_,a_] = 1/(2a)(UnitStep[x-(-a)] - UnitStep[x-a])

Plot[del[x,1],{x,-8,8}]

Integrate[del[x,1],{x,-100,100}]

We examine an animation or a Manipulation of del[x,a] as a gets smaller an smaller. We say DiracDelta[x] is the limit of del[x,a] as a->0.

Here are the plots in Manipulation to illustrate this point.  By clicking on the + box in this command one can see the values of a and by using the slider next to a one can see the effect of decreasing the value of a in the plot of the resulting function.

Manipulate[Plot[del[x,a],{x,-8,8},PlotRange->{{-8,8},{0,10}},PlotPoints->100],{a,3,.001,-.0001}]

We can combine UnitStep functions to turn things on and then off. Consider this function.

h[x_] = UnitStep[x-3] - UnitStep[x-10]

Plot[h[x],{x,0,20},PlotStyle->{Thickness[.01]}]

Suppose we have a mass spring system  and we want to see the effect of adding additional weight (force)  on the mass after it is in motion, say adding 5 units of force at time t = 2.

eq = LaplaceTransform[1 y''[t] + 6 y'[t]  + 5 y[t]==5 UnitStep[t-2] ,t,s]/.{LaplaceTransform[y[t],t,s]->Y[s],y[0]->4,y'[0]->0}

Now we solve for Y[s], the Laplace Transform of y[t] in order to get ready for the inverse Laplace Transform.

L[s_] =Y[s]/.Solve[eq,Y[s]][[1]]

Now we apply the Inverse Laplace Transform to get the solution for y[t].

yl[t_] = InverseLaplaceTransform[L[s],s,t]

From the plot we see the change in the action at t = 2 due to our applying additional force at that time.

Plot[yl[t],{t,0,10},PlotRange->{{0,10},{0,4}}]

Suppose we have a spring mass system  and we want to see the effect of hitting it with  spike force of 4 units of force after two (t = 2) seconds of "usual" motion.

eq = LaplaceTransform[1 y''[t] + 6 y'[t]  + 5 y[t]==4 DiracDelta[t-2] ,t,s]/.{LaplaceTransform[y[t],t,s]->Y[s],y[0]->0,y'[0]->0}

Now we solve for Y[s], the Laplace Transform of y[t] in order to get ready for the inverse Laplace Transform.

L[s_] =Y[s]/.Solve[eq,Y[s]][[1]]

Now we apply the transform to get the solution for y[t].

yl[t_] = InverseLaplaceTransform[L[s],s,t]

From the plot we see the change in the action at t = 2 due to our applying a spike force at that time. Up until the spring mass was at rest due to it initial conditions of y(0) = 0 and y’(0) = 0.

Plot[yl[t],{t,0,10},PlotRange->{{0,10},{0,1}}]

Let us consider a tank system in which we "dump" stuff (DiracDelta function) and we alter amounts suddenly (UnitStep function).  First the system without any drastic dumping. Can you draw a scenario for this system given the equations?

eq={x'[t]==(4*y[t])/200-(7*x[t])/200,y'[t]==(7*x[t])/200-(7*y[t])/200,x[0]==50,y[0]==0}

sol=DSolve[eq,{x[t],y[t]},t]

xs[t_]=x[t]/.sol[[1]];ys[t_]=y[t]/.sol[[1]];

Plot[{xs[t],ys[t]},{t,0,100},PlotRange->{{0,100},{0,50}}]

Let us consider a tank system in which we (1) "dump" stuff (DiracDelta function), (2) we alter amounts suddenly (UnitStep) function, or combinations of these. What are we dumping into the system?  Where are we dumping it in?  When do we start dumping?  How often are we dumping it in?

What does these stuff functions do?

stuff1[t_] = 40DiracDelta[t-25];

stuff2[t_] = 40UnitStep[t-25];

stuff3[t_] = 40UnitStep[t-25] - 20 UnitStep[t-30];

stuff4[t_] = Sum[20DiracDelta[t-a],{a,20,200,20}];

We might then modify the model system with the new "stuff" function and Laplace away!

stuff[t_] = stuff1[t]

eq = {x'[t] == 4 y[t]/200 - 7 x[t]/200 + stuff[t],
	y'[t] == 7 x[t]/200 - 7 y[t]/200,
	x[0]== 50, y[0]==0}

eq[[1]]

Now use LaplaceTransform syntax, solve, plot, and examine solution.

eq1 = LaplaceTransform[eq[[1]] ,t,s]/.{LaplaceTransform[x[t],t,s]->X[s],LaplaceTransform[y[t],t,s]->Y[s],
x[0]->50,y[0]->0}

eq2 = LaplaceTransform[eq[[2]] ,t,s]/.{LaplaceTransform[x[t],t,s]->X[s],LaplaceTransform[y[t],t,s]->Y[s],
x[0]->50,y[0]->0}

XYsol = Solve[{eq1,eq2},{X[s],Y[s]}]

Xsol[s_] = X[s]/.XYsol[[1]]//Expand

Ysol[s_] = Y[s]/.XYsol[[1]]//Expand

xsol[t_] = InverseLaplaceTransform[Xsol[s],s,t]

ysol[t_] = InverseLaplaceTransform[Ysol[s],s,t]

After solving using Laplace Transform techniques we plot our solution for x(t). We see that from the dumping due to stuff[t] = 40*DiracDelta[t-25] of 40 additional units instantly at time t = 25 there is a jump in x(t) at time t = 25 and an attendant rapid change in the second tank, y(t).

xPlot = Plot[xsol[t],{t,0,200},PlotRange->{{0,200},{0,100}}]

yPlot=Plot[ysol[t],{t,0,200},PlotRange->{{0,200},{0,100}},PlotStyle->Hue[.3]]

Here we see both x (t) in gold and y (t) in green showing the effect of our dump in both tanks.

Show[xPlot,yPlot]

Suppose we have a spring mass system  and we want to see the effect of turning on a force and turning off that force, over and over again, perhaps by placing the mass in an alternating magnetic field.   We model this  with small rectangular functions of forces, 0.5 units high and 1 unit long, spaced 2 units apart, centered at 1, 3, 5, etc.

rut[t_] = Sum[5  UnitStep[t-2 k] - 5 UnitStep[t - (2k+1)],{k,0,10}]

pRut = Plot[rut[t],{t,0,20},PlotRange->{{0,20},{0,6}}]

We apply our force function rut(t).

eq = LaplaceTransform[1 y''[t] + 6 y'[t]  + 5 y[t]==rut[t] ,t,s]/.{LaplaceTransform[y[t],t,s]->Y[s],y[0]->0,y'[0]->0}

Now we solve for Y[s], the Laplace Transform of y[t] in order to get ready for the inverse Laplace Transform.

L[s_] =Y[s]/.Solve[eq,Y[s]][[1]]

We finally apply the Inverse Laplace Transform to get the solution for y[t].

yl[t_] = InverseLaplaceTransform[L[s],s,t]

Here is our solution plot. Notice the effects of turning on and turning off the driver force.

pSol = Plot[yl[t],{t,0,20},PlotRange->{{0,20},{0,2}}]

We examine the input function (square wave) and the output function (rounded wave).

Show[pRut,pSol]

Up close. Notice the delay of the system (system response) in responding to the turn on and to the shut off of the input function

Show[pRut,pSol,PlotRange->{{0,5.5},{0,6}}]

There are many uses and reasons for using Laplace Transforms, indeed, often engineers live in the frequency domain of the transform world, i.e. s, without returning to the time domain world of t, for they learn to "see" and understand in the frequency world and, as we saw, it is algebra and not calculus - thus making it more “elementary.”