Residue Method

Now we define the inverse Laplace transform using the 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(λ), where the degree of Q(λ) must be greater than that of P(λ). 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 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 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 m_i. 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 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) . \]

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 in the following example.

Example: Let us consider 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]= -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])

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] . \]

■