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

Preface


This section presents the residue method and its application to determination of the inverse Laplace transforms. Although the residue technique is widely used in analysis and theory of functions of a complex variable, its presentation is based on the novel approach developed by Vladimir Dobrushkin. It embraces previously discussed methods: partial fraction decomposition and the convolution formula.

Residue Method


Now we define the inverse Laplace transform using the residue method, assuming that singular points are at most of multiplicity 3. It is applied to irreducible expressions of the form: F(λ) = P(λ)/Q(λ), where the degree of Q(λ) must be greater than that of P(λ). Zeroes of the denominator polynomial are called poles of the meromorphic function F(λ). Then the inverse Laplace transform of the meromorphic function F(λ) is the sum of all residues of \( e^{\lambda t} F(\lambda ) \) evaluated at the poles of F(λ).

\begin{equation} \label{EqResidue.1} 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} H(t) , \end{equation}
where the 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 a negative argument. The residues at each null are evaluated according to the following rule.

Suppose that the meromorphic 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 null 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

\begin{equation} \label{EqResidue.2} \mbox{Res}_{\lambda = \lambda_i} \, F(\lambda ) \, e^{\lambda \, t} = \frac{P(\lambda_i )}{Q' (\lambda_i )} \, e^{\lambda_i t} . \end{equation}

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

\begin{equation} \label{EqResidue.3} \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] . \end{equation}
In general, when λi is an n-fold root of \( Q(\lambda ) =0 , \) then
\begin{equation} \label{EqResidue.4} \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] . \end{equation}

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
\begin{equation} \label{EqResidue.5} \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\} , \end{equation}
where
\begin{equation} \label{EqResidue.6} 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 )} . \end{equation}

Example 1: 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 2: 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 according to Eq.\eqref{EqResidue.4}, 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*}
D[(s^2 + 4*s + 1)/(s - 3)*Exp[s*t], s, s]
(2 E^(s t))/(-3 + s) - (2 E^(s t) (4 + 2 s))/(-3 + s)^2 + ( 2 E^(s t) (1 + 4 s + s^2))/(-3 + s)^3 + ( 2 E^(s t) (4 + 2 s) t)/(-3 + s) - ( 2 E^(s t) (1 + 4 s + s^2) t)/(-3 + s)^2 + ( E^(s t) (1 + 4 s + s^2) t^2)/(-3 + s)
% /. s -> 2
-44 E^(2 t) - 42 E^(2 t) t - 13 E^(2 t) t^2
and
(s^2 + 4*s + 1)/(s - 2)^3*Exp[s*t] /. s -> 3
22 E^(3 t)
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) . \]

Finally, we check the answer with Mathematica:

InverseLaplaceTransform[(s^2 + 4*s + 1)/(s - 2)^3/(s - 3), s, t]
1/2 E^(2 t) (-44 + 44 E^t - 42 t - 13 t^2)
and
Simplify[LaplaceTransform[1/2 E^(2 t) (-44 + 44 E^t - 42 t - 13 t^2),t,s]
(1 + 4 s + s^2)/((-2 + s)^3 (-3 + s))
End of Example 2
   ■

We present Mathematica code to determine the inverse Laplace transform based on the residue method.

(* 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 3: We demonstrate an application of 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) . \]

As another option, we use the subroutine described previously:
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}}

Finally, we calculate residues manually with the aid of Mathematica.

For residue at λ = j, we have

\[ \mbox{Res}_{\lambda = {\bf j}} \frac{2\lambda -1}{(\lambda^2 +1)(\lambda -2)^3} = \lim_{\lambda \to {\bf j}} \frac{2\lambda -1}{2\lambda\,(\lambda -2)^3} = \left( \frac{7}{250} - \frac{12{\bf j}}{125} \right) e^{{\bf j}t} . \]
Assuming[ t > 0, Exp[s*t]*(2*s - 1)/2/s/(s - 2)^3 /. s -> I]
(7/250 - (12 I)/125) E^(I t)
We don't need to evaluate the residue at λ = −j because its value is a complex conjugate of the residure at λ = j. So their sum is just real part times 2:
\[ \left( \mbox{Res}_{\lambda = {\bf j}} + \mbox{Res}_{\lambda = -{\bf j}} \right) \frac{2\lambda -1}{(\lambda^2 +1)(\lambda -2)^3} \,e^{\lambda t} = \frac{7}{125}\,\cos t + \frac{24}{125}\,\sin t . \]
We calculate the residure at the tripple pole λ = 2 using Eq.\eqref{EqResidue.4}
\[ \mbox{Res}_{\lambda = 2} \frac{(2\lambda -1)\,e^{\lambda t}}{(\lambda^2 +1)(\lambda -2)^3} = \frac{1}{2}\,\lim_{\lambda \to 2} \frac{{\text d}^2}{{\text d}\lambda^2} \, \frac{2\lambda -1}{(\lambda^2 +1)} \,e^{\lambda t} = e^{2t} \left( \frac{3}{5}\, t^2 - \frac{4}{25}\, t - \frac{14}{125} \right) . \]
D[(2*s - 1)/(s^2 + 1)*Exp[s*t], s, s]
(8 E^(s t) s^2 (-1 + 2 s))/(1 + s^2)^3 - (8 E^(s t) s)/(1 + s^2)^2 - ( 2 E^(s t) (-1 + 2 s))/(1 + s^2)^2 - ( 4 E^(s t) s (-1 + 2 s) t)/(1 + s^2)^2 + (4 E^(s t) t)/(1 + s^2) + ( E^(s t) (-1 + 2 s) t^2)/(1 + s^2)
% /. s -> 2
-((14 E^(2 t))/125) - 4/25 E^(2 t) t + 3/5 E^(2 t) t^2
   ■
There is another option, which we demonstrate in the following example.

Example 4: Let us consider the function:

\[ F(\lambda ) = \frac{2\lambda -1}{(\lambda^2 +1)^2 (\lambda -2)^3} . \tag{4.1} \]
Now, we have two double complex conjugate roots \( \lambda = \pm {\bf j} . \) The denominator Q(λ) = (λ² + 1)²(λ - 2)³ also has a triple real root λ = 2. The first step is to define the denominator Q(λ) and the numerator P(λ) = 2λ −1 in Mathematica:

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
   ■

 

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)