Laplace Transformation

This tutorial contains many Mathematica scripts. You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to <Vladimir_Dobrushkin@brown.edu>

Return to computing page
Return to computing page for the second course APMA0330
Return to Mathematica tutorial for the first course
Return to Mathematica tutorial for the second course
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340

Contents

Part VI: Laplace Transformation

Laplace transform
Heaviside function
Laplace Transform of Discontinuous Functions
Inverse Laplace transformation
Laplace transformation in differential equations

Mechanical and Electrical Vibrations
Other applications

1.6.1. Laplace Tranform

LaplaceTransform[t^4 Sin[t], t, lambda]
LaplaceTransform[HeavisideTheta[t - 1] t, t, lambda]
LaplaceTransform[DiracDelta[t], t, lambda]
LaplaceTransform[f''[t], t, lambda]

Out[4]= -lambda f[0] + lambda^2 LaplaceTransform[f[t], t, lambda] - Derivative[1][f][0]

LaplaceTransform[If[t > Pi, Cos[t], t], t, s]

Out[4]= (E^(-\[Pi] s) (-1 + E^(\[Pi] s) - \[Pi] s - s^3/(1 + s^2)))/s^2

LaplaceTransform[UnitStep[2 Sin[t]], t, lambda]

Out[5]= E^(\[Pi] lambda)/((1 + E^(\[Pi] lambda)) lambda)

LaplaceTransform[ArcTan[t/a], t, s, Assumptions -> a > 0]

Out[5]= (2 CosIntegral[a s] Sin[a s] + Cos[a s] (\[Pi] - 2 SinIntegral[a s]))/(2 s)

LaplaceTransform[y''[t] + 4 y[t] == 0, t, lambda]

Out[4]= 4 LaplaceTransform[y[t], t, lambda] + lambda^2 LaplaceTransform[y[t], t, lambda] - lambda y[0] - Derivative[1][y][0] == 0

% /. {y[0] -> 1, y'[0] -> 2}

Out[5]= -2 - lambda + 4 LaplaceTransform[y[t], t, lambda] + lambda^2 LaplaceTransform[y[t], t, lambda] == 0

Solve[%, LaplaceTransform[y[t], t, lambda]]

Out[6]= {{LaplaceTransform[y[t], t, lambda] -> (2 + lambda)/(4 + lambda^2)}}

InverseLaplaceTransform[%, lambda, t]

Out[7]= {{y[t] -> Cos[2 t] + Sin[2 t]}}

Now we give an example of application of the inverse Laplace transform.

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

Out[8]= 4 (1/4 - 1/4 Cos[2 t])

Simplify[%]

Out[9]= 2 Sin[t]^2

Apart[4/s/(4 + s^2), s]

Out[10]= 1/s - s/(4 + s^2)

LaplaceTransform[2 Sin[t]^2, t, lambda]

Out[11]= 4/(lambda (4 + lambda^2))

1.6.2. Heaviside Function

I. How to define Heaviside functions

To define a function, just type in the formula. We need to use a special form for the left hand side, which includes an underscore after the name of the variable, and a special "colon-equals" sign for the function definition:

f[x_] := (Cos[x] -1)/ x^2

There is no output on this imput. To see it, type Print[f[x]] You can use this function with different arguments or obtain its numerical values:

f[2 x+1]

Out[2]= (Cos[2 x+1]-1)/(2 x+1)^2

N[f[1.5]]

Out[3]= -0.41306

The simplest user-defined functions are the "one-liners", where the quantity of interest can be computed by a single formula. However, in may cases, you may find it impossible to define the function's value in a single simple formula. Instead, you may need to carry out several steps of computation, using temporary variables. You may want several input values, and you may want the user to group some of those input values in curly brackets.

1.6.3. Laplace Transform of Discontinuous Functions

I. How to define functions

g[t_] = Piecewise[{{t, 0 < t < 1}, {(t - 2)^2, 1 < t < 2}, {6 - t^2, 2 < t}}]

LaplaceTransform[g[t], t, lambda]

Out[14]= (1/(lambda^3))E^(-2 lambda) (-4 - 4 lambda - E^lambda lambda + E^(2 lambda) lambda + 2 lambda^2 - E^lambda lambda^2 + Gamma[3, -lambda])

Shift Rule:

LaplaceTransform[UnitStep[t - 2]*Sin[t], t, lambda]

Out[5]= (E^(-2 lambda) (Cos[2] + lambda Sin[2]))/(1 + lambda^2)

LaplaceTransform[UnitStep[t - 2]*Sin[t - 2], t, lambda]

Out[6]= E^(-2 lambda)/(1 + lambda^2)

LaplaceTransform[DiracDelta[t - 2]*Sin[t - 2], t, lambda]

Out[7]= 0

LaplaceTransform[DiracDelta[t - 2]*Sin[t], t, lambda]

Out[8]= E^(-2 lambda) Sin[2]

1.6.4. Inverse Laplace Transform

I. Standard procedure

One of the methods to find inverse laplace transform is to us epartial fraction decomposition and then apply inverse Laplace transform to each term (usually using a table of given Laplace transforms). However, Mathematica allows one to find the inverse Laplace transform stright forward:

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

Out[8]= 4 (1/4 - 1/4 Cos[2 t])

Simplify[%]

Out[9]= 2 Sin[t]^2

Apart[4/s/(4 + s^2), s]

Out[10]= 1/s - s/(4 + s^2)

LaplaceTransform[2 Sin[t]^2, t, lambda]

Out[11]= 4/(lambda (4 + lambda^2))

The general approach:

InverseLaplaceTransform[%, lambda, t]

Out[7]= {{y[t] -> Cos[2 t] + Sin[2 t]}}

II. Convolution rule

III. Residue method

Now we define the inverse Laplace transform using residue method, assuming that singular points are at most of multiplicity 3. It is applyied to the expression of the form: p(s)/q(s). The output is a collection of residues for each singular point.

(* use variable "s" to define p and q *)

Natalia[p_, q_] :=
(g = ToRules[Roots[q == 0, s]];
u = s /. {g};
h = Tally[u];
j = Length[r];
Res1 = {s - s0}*{p/{Factor[q, Extension -> I]}}*Exp[s*t];
Res2[s_] = {s - s0}^2*{p/{Factor[q, Extension -> I]}}*Exp[s*t];
Res3[s_] = {s - s0}^3*{p/{Factor[q, Extension -> I]}}*Exp[s*t];
r = Table[h[[i, 2]], {i, 1, Length[h]}];
Do[If[r[[i]] == 1, Res = Res1 /. s0 -> h[[i, 1]];
Res = ComplexExpand[Res /. s -> h[[i, 1]]]; Print[Res],
If[r[[i]] == 2, ResII[s_] = Res2[s] /. s0 -> h[[i, 1]];
ResII[s_] = ComplexExpand[ResII'[s] /. s -> h[[i, 1]]];
Print[ResII[s]],
If[r[[i]] == 3, ResIII[s_] = Res3[s] /. s0 -> h[[i, 1]];
ResIII[s_] = {1/2}*{ComplexExpand[
ResIII''[s] /. s -> h[[i, 1]]]}; Print[ResIII[s]]]]], {i, 1,
j}])

Example.

Natalia[2*s - 1, (s^2 + 1)*(s - 2)^3]

{{{1/2 (-((14 E^(2 t))/125)-4/25 E^(2 t) t+3/5 E^(2 t) t^2)}}}
{{(7 Cos[t])/250+I (-((12 Cos[t])/125)+(7 Sin[t])/250)+(12 Sin[t])/125}}
{{(7 Cos[t])/250+I ((12 Cos[t])/125-(7 Sin[t])/250)+(12 Sin[t])/125}}

There is another option, which we demonstrate on the same example.

P = 2 s - 1 (* Defines polynomial P *)

Out[1]= -1 + 2 s

Q = (s^2 + 1)^2 (s - 2)^3 (* Defines polynomial Q *)

Out[2]= (-2 + s)^3 (1 + s^2)^2

s /. Solve[Q == 0] (* Finds the roots of Q(s)=0 *)

Out[3]= {-I, -I, I, I, 2, 2, 2}

Tally[%] (* Finds the number of repeats of each root (multiplicity)*)

Out[4]= {{-I, 2}, {I, 2}, {2, 3}}

MatrixForm[%] (* Creates solutions matrix with roots in one column and multiplicities in another *)

Out[5]=
(-i 2 )
( i 2 )
( 2 3 )

Dimensions[%] (* Measures the dimensions of the matrix (2x2),(3x2),(6x2), etc*)

Out[6]= {3, 2}

v = Part[%, 1] (* Finds the length of the matrix (v=number of rows) *)

Out[7]= 3

Here we start the Do loop to find residue for each root, where each row of the solutions matrix contains a root and its multiplicity

Do[T := Tally[s /. Solve[Q == 0]];
M := MatrixForm[T];    (* Recreates the solutions matrix for use in Do loop *)
r := Extract[Extract[Extract[M, 1], l], 1];  (* Extracts the root (r) from row of the matrix *)
m := Extract[Extract[Extract[M, 1], l], 2]; (* Extracts multiplicity (m) of root from row of the matrix *)
f = ((((s - r)^m) (E^(s*t)) (P))/Q);
F = Simplify[f];    (* Simplifies f into nicely factored form *)
d = (1/(m - 1)!) D[F, {s, If[m > 1, m - 1, 0]}];  (* it is the equation for residue *)
Res = d /. s -> r;  (* Calculates residue at a particular root (r=root) *)
Print[Y = ComplexExpand[Re[Res], {I}]], {l, v}]

It takes the Euler expansion of previous line of code and extracts all real parts of that expression. Final output will contain several terms, that are functions of t, that will be the total residue when added. Note that the Res function computes the residue for each root by taking roots and multiplicities from row, where l ranges from 1 to v, the length of the solutions matrix.

Out[8]= -((17 Cos[t])/625)-6/125 t Cos[t]+(31 Sin[t])/625+7/500 t Sin[t] -((17 Cos[t])/625)-6/125 t Cos[t]+(31 Sin[t])/625+7/500 t Sin[t] (34 E^(2 t))/625-14/125 E^(2 t) t+3/50 E^(2 t) t^2

To check the answer, we type:

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

Out[9]=
(1/1250)(68 E^(2 t) - 140 E^(2 t) t + 75 E^(2 t) t^2 - 68 Cos[t] -
120 t Cos[t] + 124 Sin[t] + 35 t Sin[t])

MATHEMATICA TUTORIAL for the First Cource. Part VI.

Ver. 1.0

Vladimir Dobrushkin

Laplace Transformation

Contents [hide]

1.6.1. Laplace Tranform

1.6.2. Heaviside Function

I. How to define Heaviside functions

1.6.3. Laplace Transform of Discontinuous Functions

I. How to define functions

1.6.4. Inverse Laplace Transform

I. Standard procedure

II. Convolution rule

III. Residue method

1.6.5. Differential Equations

I. How to define functions

I.Differential Equations with Discontinuous Functions

1.6.6. Mechanical and Electrical Applications

6.7. References