# Preface

This tutorial was made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

# Residue Method

Now we define the inverse Laplace transform using residue method, assuming that singular points are at most of multiplicity 3. Although the residue method is widely used in complex analysis, it was adapted by Vladimir Dobrushkin to embrace previously discussed methods: partial fraction decomposition and the convolution formula. It is applied to the expression of the form: P(λ)/Q(λ). The output is a collection of residues for each singular point.

$f(t) = {\cal L}^{-1} \left[ F(\lambda ) \right] = {\cal L}^{-1} \left[ \frac{P(\lambda )}{Q(\lambda )} \right] = \sum_{i=1}^N \,\mbox{Res}_{\lambda = \lambda_i} \, F(\lambda ) \, e^{\lambda \, t} ,$
where sum ranges over all zeroes of the equation $$Q(\lambda ) =0 ,$$ and the function f must be multiplied by the Heaviside function to ensure that f(t) is identically zero for negative argument. The residures at each null are evaluated according to the following rule.

Suppose that a function $$F(\lambda ) = \frac{P(\lambda )}{Q(\lambda )}$$ is a ratio of two irreducible polynomials (or entire functions). Let $$\lambda = \lambda_i , \ i=1,2,\ldots , N$$ be nulls of the denominator ($$N = \infty$$ if Q(λ) is an entire function). Recall that if $$Q(\lambda_i ) =0, Q' (\lambda_i ) =0, \ldots , Q^{(m_i -1)} (\lambda_i ) =0, Q^{(m_i )} (\lambda_i \ne 0 ,$$ then we say that the denominator has zero at $$\lambda = \lambda_i$$ of multiplicity mi. In this case, we say that the function F(λ) has poles at $$\lambda = \lambda_i , \quad i=1,2,\ldots , N .$$ Note that the numerator is not zero at these poles.

If $$\lambda = \lambda_i$$ is a simple root (meaning that its multiplicity is 1) of the equation $$Q(\lambda ) =0 ,$$ then

$\mbox{Res}_{\lambda = \lambda_i} \, F(\lambda ) \, e^{\lambda \, t} = \frac{P(\lambda_i )}{Q' (\lambda_i )} \, e^{\lambda_i t} .$

If $$\lambda = \lambda_i$$ is a double root (meaning that $$Q(\lambda_i ) = Q' (\lambda_i ) =0, \quad\mbox{but} \quad Q'' (\lambda_i ) \ne 0 ,$$ ), then

$\mbox{Res}_{\lambda = \lambda_i} \, F(\lambda ) \, e^{\lambda \, t} = \lim_{\lambda \mapsto \lambda_i} \, \frac{\text d}{{\text d}\lambda} \left[ (\lambda - \lambda_i )^2 \,\frac{P(\lambda_i )}{Q (\lambda_i )} \, e^{\lambda_i t} \right] .$
In general, when λi is an n-fold root of $$Q(\lambda ) =0 ,$$ then
$\mbox{Res}_{\lambda = \lambda_i} \, F(\lambda ) \, e^{\lambda \, t} = \frac{1}{(n-1)!} \,\lim_{\lambda \mapsto \lambda_i} \, \frac{{\text d}^{n-1}}{{\text d}\lambda^{n-1}} \left[ (\lambda - \lambda_i )^n \,\frac{P(\lambda_i )}{Q (\lambda_i )} \, e^{\lambda_i t} \right] .$

Note that if $$\lambda_i = \alpha + {\bf j} \beta$$ is a complex null of the denominator Q(λ) (the root of the polynomial equation with real coefficients $$Q(\lambda ) =0$$ ), then $$\overline{\lambda_i} = \alpha - {\bf j} \beta$$ is also a null of Q(λ). In this case, we don't need to calculate the residue at the complex conjugate because we know that the output will be a complex conjugate of the residue at $$\alpha + {\bf j} \beta .$$ Therefore, if the denominator Q(λ) has two complex conjugate roots $$\lambda = \alpha \pm {\bf j} \beta ,$$ then the sum of residues at these points is just double value of the real part of one of them:

$\mbox{Res}_{\lambda = \alpha + {\bf j} \beta} \, F(\lambda ) \, e^{\lambda \, t} + \mbox{Res}_{\lambda = \alpha - {\bf j} \beta} \, F(\lambda ) \, e^{\lambda \, t} = 2\,\Re \, \mbox{Res}_{\lambda = \alpha + {\bf j} \beta} \, F(\lambda ) \, e^{\lambda \, t} .$
In the case of simple complex roots, we get
$\mbox{Res}_{\lambda = \alpha + {\bf j} \beta} \, F(\lambda ) \, e^{\lambda \, t} + \mbox{Res}_{\lambda = \alpha - {\bf j} \beta} \, F(\lambda ) \, e^{\lambda \, t} = 2\,e^{\alpha \,t} \left\{ A\,\cos \beta t - B\,\sin \beta t \right\} ,$
where
$A = \Re \,\mbox{Res}_{\lambda = \alpha + {\bf j} \beta} \, F(\lambda ) = \mbox{Re} \, \frac{P(\alpha + {\bf j} \beta )}{Q' (\alpha + {\bf j} \beta )} , \qquad B = \Im \,\mbox{Res}_{\lambda = \alpha + {\bf j} \beta} \, F(\lambda ) = \mbox{Im} \, \frac{P(\alpha + {\bf j} \beta )}{Q' (\alpha + {\bf j} \beta )} .$

Example: Let us find the inverse Laplace transform of the function

$F(\lambda ) = \frac{\lambda^2 +4\,\lambda +1}{(\lambda +1)(\lambda -2)(\lambda -3)}$
The denominator has three simple roots $$\lambda = -1, \ 2, \ 3 .$$ Therefore, its inverse Laplace transform is the sum of residues:
${\cal L}^{-1} \left[ \frac{\lambda^2 +4\,\lambda +1}{(\lambda +1)(\lambda -2)(\lambda -3)} \right] = \left( \mbox{Res}_{\lambda =-1} + \mbox{Res}_{\lambda = 2} + \mbox{Res}_{\lambda =3} \right) F(\lambda )\, e^{\lambda \,t} .$
We calculate each residue separately.
\begin{align*} \mbox{Res}_{\lambda =-1} \, F(\lambda )\, e^{\lambda \,t} &= \lim{\lambda \mapsto -1} \,\frac{\lambda^2 +4\,\lambda +1}{(\lambda -2)(\lambda -3)} \, e^{\lambda\,t} = -\frac{1}{6}\, e^{-t} , \\ \mbox{Res}_{\lambda =2} \, F(\lambda )\, e^{\lambda \,t} &= \lim{\lambda \mapsto 2} \,\frac{\lambda^2 +4\,\lambda +1}{(\lambda +1)(\lambda -3)} = -\frac{13}{3}\, e^{2t} , \\ \mbox{Res}_{\lambda =3} \, F(\lambda )\, e^{\lambda \,t} &= \lim{\lambda \mapsto 3} \,\frac{\lambda^2 +4\,\lambda +1}{(\lambda +1)(\lambda -2)} = \frac{11}{2}\, e^{3t} . \end{align*}
Adding all residues, we find the inverse Laplace transform
${\cal L}^{-1} \left[ \frac{\lambda^2 +4\,\lambda +1}{(\lambda +1)(\lambda -2)(\lambda -3)} \right] = \left( \frac{11}{2}\, e^{3t} - \frac{13}{3}\, e^{2t} - \frac{1}{6}\, e^{-t} \right) H(t) .$

Example: Let us find the inverse Laplace transform of the function that has multiple poles

$F(\lambda ) = \frac{\lambda^2 +4\,\lambda +1}{(\lambda -2)^3 (\lambda -3)} ,$
where singular point $$\lambda =2$$ has multiplicity 3. Calculating residues, we get
\begin{align*} \mbox{Res}_{\lambda =2} \,F(\lambda ) \, e^{\lambda\, t} & = \frac{1}{2}\,\lim_{\lambda \mapsto 2} \, \frac{{\text d}^2}{{\text d} \lambda^2} \, \frac{\lambda^2 +4\,\lambda +1}{(\lambda -3)} \, e^{\lambda \,t} = - \left( 22+ 21\,t + \frac{13}{2}\,t^2 \right) e^{2t} , \\ \mbox{Res}_{\lambda =3} \,F(\lambda ) \, e^{\lambda\, t} & = \lim_{\lambda \mapsto 3} \, \frac{\lambda^2 +4\,\lambda +1}{(\lambda -2)^3} \, e^{\lambda\, t} = 22\, e^{3t} . \end{align*}
Adding these residues, we find the inverse Laplace transform
${\cal L}^{-1} \left[ \frac{\lambda^2 +4\,\lambda +1}{(\lambda -2)^3 (\lambda -3)} \right] = 22\, e^{3t} \, H(t) - e^{2t} \left( 22+ 21\,t + \frac{13}{2}\,t^2 \right) H(t) .$

We present Mathematica code to calculate the residue.

(* 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: We demonstrate application of the our code to find the inverse Laplace transform of

$F(\lambda ) = \frac{2\lambda -1}{(\lambda^2 +1)(\lambda -2)^3} .$
First, we use the standard Mathematica command
InverseLaplaceTransform[2*s-1)/(s^2+1)/(s-2)^3 , s, t]
to get
${\cal L}^{-1} \left[ \frac{2\lambda -1}{(\lambda^2 +1)(\lambda -2)^3} \right] = \frac{1}{125} \left[ 7\,\cos t + 24\,\sin t \right] H(t) - \frac{1}{250} \, e^{2t} \left[ 14 + 20t +75t^2 \right] H(t) .$
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 a similar example for the function

$F(\lambda ) = \frac{2\lambda -1}{(\lambda^2 +1)^2 (\lambda -2)^3} .$
Now we have two double complex conjugate roots $$\lambda = \pm {\bf j} .$$

P = 2 s - 1   (* Defines polynomial P *)
Out= -1 + 2 s
Q = (s^2 + 1)^2 (s - 2)^3       (* Defines polynomial Q *)
Out= (-2 + s)^3 (1 + s^2)^2
s /. Solve[Q == 0]  (* Finds the roots of Q(s)=0 *)
Out= {-I, -I, I, I, 2, 2, 2}
Tally[%]     (* Finds the number of repeats of each root (multiplicity)*)
Out= {{-I, 2}, {I, 2}, {2, 3}}
MatrixForm[%]  (* Creates solutions matrix with roots in one column and multiplicities in another *)
Out=
(-i 2 )
( i 2 )
( 2 3 )
Dimensions[%]   (* Measures the dimensions of the matrix (2x2),(3x2),(6x2), etc*)
Out= {3, 2}
v = Part[%, 1]   (* Finds the length of the matrix (v=number of rows) *)
Out= 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= -((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=
(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])

To check, we make manual calculations.
\begin{align*} \mbox{Res}_{\lambda = 2} \, F(\lambda )\, e^{\lambda\,t} &= \left. \frac{1}{2}\,\frac{{\text d}^2}{{\text d}\lambda^2} \, \frac{2\lambda -1}{(\lambda^2 +1)}\, e^{\lambda\,t} \right\vert_{\lambda =2} = \lim_{\lambda \mapsto 2} \frac{e^{\lambda\,t}}{2(1+\lambda^2 )^3} \left[ 2 + 4t-t^2 + 2\lambda^5 t^2 -\lambda^4 t(4+t) -2\lambda^2 (3+t^2) +4\lambda^3 (1+t+t^2 ) +2\lambda (2t+t^2 -6) \right] \\ &= \frac{1}{250}\, e^{2t} \left[ 75t^2 -20t -14 \right] , \\ \mbox{Res}_{\lambda = {\bf j}} \, F(\lambda )\, e^{\lambda\,t} &= \left. \frac{{\text d}}{{\text d}\lambda} \,\frac{2\lambda -1}{(\lambda -2)^3}\, e^{\lambda\,t} \right\vert_{\lambda ={\bf j}} = \lim_{\lambda \mapsto {\bf j}} \frac{e^{\lambda\,t}}{(\lambda -2)^4} \left[ -1 + 2t + 2\lambda^2 t - \lambda (4+5t) \right] \\ &= \frac{24{\bf j} -7}{625} \, e^{{\bf j}t} \left[ -1 - {\bf j} (4+5t) \right] . \end{align*}
Extracting the real part and multiplying by 2, we obtain
$2\,\mbox{Re}\, \mbox{Res}_{\lambda = {\bf j}} \, F(\lambda )\, e^{\lambda\,t} = \frac{1}{625} \left[ 14\,\cos t -14 (4+5t) \,\sin t + 48\,\sin t + 48\,(4+5t)\,\cos t \right] .$