On the Approximate Inversion of the Laplace Transform

V. M. Osipov

(Presented by the scientific-methodological seminar of the Department of Engineering and Computational Mathematics).

In many cases, inversion of the Laplace transform by the classical scheme, i.e., by computing the Bromwich contour integral, is difficult to carry out. Thus, when the Laplace transform is given in the form of a rational fraction, inversion reduces to determining the roots of an algebraic equation, which is a technically difficult task, especially for a high degree of the equation. Moreover, the problem often arises of approximating transcendental operational images by more convenient expansions in the form of rational fractions.

Below we present methods for the approximate solution of the inversion and approximation problems that are practically convenient and sufficiently effective for a broad class of operational images.

1. Inversion by means of Chebyshev “e”-functions of the third kind

The exponential Chebyshev functions of the third kind [4]

\[ S_n(t)=2e^{-\frac{at}{2}}\sqrt{1-e^{-at}}\,U^*_{n-1}(t) \qquad (n=1,2,\ldots), \tag{1-1} \]

where \(U^*_{n-1}(t)\) \((n=1,2,\ldots)\) are the Chebyshev “e”-polynomials of the second kind, being eigenfunctions of a self-adjoint differential operator, form on the interval \((0,\infty)\) a complete orthogonal system with weight factor

\[ W(t)=e^{-\frac{at}{2}}(1-e^{-at})^{-\frac12}. \]

In other words, the series

\[ f(t)=\sum_{k=1}^{\infty} B_k S_k(t), \tag{1-2} \]

where

\[ B_k=\frac{2a}{\pi}\int_0^\infty \frac{e^{-\frac{at}{2}}f(t)}{\sqrt{1-e^{-at}}}\,S_k(t)\,dt \qquad (k=1,2,\ldots) \tag{1-3} \]

converges at every interior point that is a point of continuity of \(f(t)\), if the condition is satisfied

\[ \int_0^\infty W(t)[f(t)]^2\,dt < \infty . \tag{1-4} \]

At discontinuities of the first kind \(t=t_i\) \((i=1,2,\ldots,n)\), the series (1-2) converges to the value

\[ \frac{1}{2}\,[f(t_i+0)-f(t_i-0)]. \]

The coefficients \(B_k\) are very simply connected with the Laplace transform of the function \(f(t)\) being expanded.

Let us introduce the concept of the exponential moment of the function \(f(t)\). This is a real quantity defined by the expression

\[ M^*_{k-1}=\int_0^\infty e^{-kat}f(t)\,dt \qquad (k=1,2,\ldots). \tag{1-5} \]

Let \(F(p)\) be the Laplace transform of the function \(f(t)\), i.e.,

\[ F(p)=\int_0^\infty e^{-pt}f(t)\,dt. \]

Comparing with (1-5), we arrive at the conclusion that

\[ M^*_{k-1}=F(ka) \quad \text{or} \quad M^*_k=F[(k+1)a]. \tag{1-6} \]

Thus, the exponential moment of order \(k\) of the function \(f(t)\) is equal to the Laplace transform of this function at \(p=(k+1)a\). Substituting (1-1) into (1-3), we obtain the following formula for the coefficient \(B_k\):

\[ B_k=\frac{4a}{\pi}\int_0^\infty e^{-at}f(t)U^*_{k-1}(t)\,dt. \tag{1-7} \]

The explicit expression for the “\(e\)-polynomial” \(U^*_{k-1}(t)\) has the form [4]

\[ U^*_{k-1}(t)=\sum_{n=0}^{k-1}\beta_{kn}e^{-nat}, \]

where

\[ \tag{1-8} \]

\[ \beta_{kn}=\frac{(-1)^{n+k-1}2^n(n+k)!}{(2n+1)(k-1-n)!\,n!\,(2n-1)!!}. \tag{1-9} \]

Thus,

\[ B_k=\frac{4a}{\pi}\sum_{n=0}^{k-1}\beta_{kn}M^*_n \qquad (k=1,2,3,\ldots), \tag{1-10} \]

i.e., the coefficient \(B_k\) is a linear combination of the first \(k\) exponential moments. Formula (1-10) can be written in matrix form

\[ \bar{B}=\frac{4a}{\pi}\,\beta \bar{M}^*, \tag{1-11} \]

where the vector \(B=\{B_1,B_2,\ldots,B_m,\ldots\}\), and the vector \(M^*=\{M^*_0,M^*_1,\ldots,M^*_m,\ldots\}\).

The coefficients of the triangular matrix \(\beta\) up to \(k=8\) inclusive are given in Appendix III; they are also the coefficients of the “\(e\)-polynomials” \(U^*_{k-1}(t)\) \((k=1,2,\ldots,8)\). The series (1-2) will converge considerably faster if \(f(0)=f(\infty)=0\). This condition means that the image \(F(p)\) of the function \(f(t)\) has no poles at the origin or on the imaginary axis and, moreover, there exists the limit

\[ \lim_{p\to\infty} p F(p)=0. \]

Let \(\Phi(p)\) be the Laplace transform of the function \(\varphi(t)\). In the general case \(\varphi(0)\ne 0\) and \(\varphi(\infty)\ne 0\); then, instead of (1-2), one should write

\[ \varphi(t)=\varphi(0)e^{-\frac{at}{2}}+\varphi(\infty)\left(1-e^{-\frac{at}{2}}\right)+\sum_{k=1}^{\infty} B_k S_k(t). \tag{1-12} \]

Denote

\[ f(t)=\varphi(t)-\varphi(0)e^{-\frac{at}{2}}-\varphi(\infty)\left(1-e^{-\frac{at}{2}}\right), \tag{1-13} \]

or, in operator form \([\Phi(p)\doteq \varphi(t)]\),

\[ F(p)=\Phi(p)-\frac{\varphi(0)}{p+\frac{a}{2}}-\frac{\frac{a}{2}\,\varphi(\infty)}{p\left(p+\frac{a}{2}\right)}, \tag{1-14} \]

where \(f(0)=f(\infty)=0\). Thus the matter has been reduced to zero boundary conditions.

Let us find the Laplace transform of the “e”-function \(S_k(t)\). Since the substitution

\[ e^{-\frac{at}{2}}=\cos \frac{\alpha}{2},\quad \text{i.e.}\quad \alpha=2\arccos e^{-\frac{at}{2}}, \]

transforms \(S_k(t)\) into \(\sin k\alpha\), and \(T_k^*(t)\) into an “e”-Chebyshev polynomial of the first kind in \(\cos k\alpha\) [4], we easily obtain

\[ \frac{dS_k(t)}{dt}=ka\,\frac{e^{-\frac{at}{2}}}{\sqrt{1-e^{-at}}}\,T_k^*(t). \tag{1-15} \]

We have [4]

\[ \frac{e^{-\frac{at}{2}}}{\sqrt{1-e^{-at}}}\,T_k^*(t)\doteq \frac{\sqrt{\pi}\,\Gamma\left(\frac{p}{a}+\frac{1}{2}\right)} {a\Gamma\left(\frac{p}{a}\right)} \frac{\prod\limits_{n=1}^{k-1}\left(\frac{p}{a}-n\right)} {\prod\limits_{n=1}^{k}\left(\frac{p}{a}+n\right)}. \]

Thus,

\[ S_k(t)\doteq (-1)^{k-1} \frac{ka\sqrt{\pi}\,\Gamma\left(\frac{p}{a}+\frac{1}{2}\right)} {p\Gamma\left(\frac{p}{a}\right)} \frac{\prod\limits_{n=1}^{k-1}(na-p)} {\prod\limits_{n=1}^{k}(na+p)}. \tag{1-16} \]

Here

\[ \prod_{n=1}^{0}\equiv 1. \]

Taking the Laplace transform of the series (1-2), we may write

\[ F(p)=\sum_{k=1}^{\infty} B_k S_k(t) = \frac{a\sqrt{\pi}\,\Gamma\left(\frac{p}{a}+\frac{1}{2}\right)} {p\Gamma\left(\frac{p}{a}\right)} \sum_{k=1}^{\infty}(-1)^{k-1}kB_k \times \]

\[ \times \frac{\prod\limits_{n=1}^{k-1}(na-p)} {\prod\limits_{n=1}^{k}(na+p)}. \tag{1-17} \]

Let us now put \(p=ma\) \((m=1,2,\ldots)\); then we shall have an infinite system of linear equations for determining \(B_k\)

\[ M^*_{m-1}=F(ma)= \frac{\sqrt{\pi}\,\Gamma\left(m+\frac{1}{2}\right)}{am\Gamma(m)} \sum_{k=1}^{\infty}(-1)^{k-1}kB_k \frac{\prod\limits_{n=1}^{k-1}(n-m)} {\prod\limits_{n=1}^{k}(n+m)}, \]

which, in expanded form, assumes a convenient triangular form, since
\(\prod\limits_{n=1}^{k-1}(n-m)=0\) for \(m\leq k-1\). Taking into account that

\[ \frac{\sqrt{\pi}\,\Gamma\left(m+\frac{1}{2}\right)}{am\Gamma(m)} = \frac{\pi}{4a}\, \frac{1\cdot3\cdot5\cdots(2m-1)}{2^{m-2}m!}, \]

we obtain

\[ \frac{4a}{\pi}M^*_0=B_1 \]

\[ \frac{4a}{\pi}M^*_1=\frac{1}{2}B_1+\frac{1}{4}B_2 \]

\[ \frac{4a}{\pi}M^*_2=\frac{5}{16}B_1+\frac{1}{4}B_2+\frac{1}{16}B_3 \]

\[ \frac{4a}{\pi}M^*_3=\frac{7}{32}B_1+\frac{7}{32}B_2+\frac{3}{32}B_3+\frac{1}{64}B_4 \]

or, in matrix form,

\[ \frac{4a}{\pi}\,\bar M^*=\beta^{-1}\bar B, \]

whence (1-11) follows. Thus, the coefficient matrix of the system (1-18) or (1-19) is the matrix inverse to the previously introduced matrix \(\beta\) with elements (1-9). In other words, the Fourier coefficients of the orthogonal expansion (1-2) can be determined as the result of an interpolation process in the complex domain, in which the interpolation nodes are the zeros of the operational image \(S_m(p)\), i.e., the zeros of the first of the discarded terms in the truncated operational series (1-17).

2. Inversion by means of integral Legendre “e”-polynomials

Inversion in the form of a segment of a series in Chebyshev exponential functions of the third kind contains an irrational factor of the form \(\sqrt{1-e^{at}x}\), which complicates analytic operations with this expression. The use, for purposes of inversion, of Legendre “e”-polynomials \(P_n^*(t)\) [4] proves more convenient from this point of view. However, the expansion in the “e”-polynomials \(P_n^*(t)\) converges considerably more slowly than the expansion (1-2). This drawback can be appreciably weakened if one uses the so-called integral Legendre “e”-polynomials \(V_n^*(t)\), which make it possible to reconstruct the function much more accurately from a given number of exponential moments, since the boundary conditions are satisfied independently of the expansion coefficients.

The integral Legendre polynomials considered by us in [3],

defined on the interval \((0,1)\); however, by the simple substitution \(x=e^{-at}\) they are transformed into “e”-polynomials \(V_n^*(t)\)

\[ V_n^*(x)\to V_n^*(e^{-at})=V_n^*(t),\qquad 0\leq t\leq \infty . \]

The basic functional relations retain their form; however, it is necessary to take into account that

\[ \frac{dV_n^*(t)}{dt}=\frac{dV_n^*(x)}{dx}\cdot \left.\frac{dx}{dt}\right|_{x=e^{-at}} . \]

In particular, we shall have [3]

\[ \left. \begin{aligned} \frac{dV_n^*(t)}{dt}&=-ae^{-at}P_n^*(t),\\ V_n^*(t)&=-a\int_0^t e^{-at}P_n^*(t)\,dt \end{aligned} \right\}, \tag{2-1} \]

\[ 2(2n+1)V_n^*(t)=P_{n-1}^*(t)-P_{n+1}^*(t), \tag{2-2} \]

\[ V_n^*(t)=e^{-at}(1-e^{-at})F(-n+1,\ n+2;\ 2;\ e^{-at}), \tag{2-3} \]

\[ V_n^*(t)=\frac{1-e^{-at}}{an(n+1)}\,\frac{dP_n^*(t)}{dt}, \tag{2-4} \]

where \(P_n^*(t)\) is the “e”-Legendre polynomial [4].

The polynomials \(V_n^*(t)\) satisfy the following differential equation:

\[ (1-e^{-at})\left[\frac{d^2V_n^*(t)}{dt^2}+a\frac{dV_n^*(t)}{dt}\right] +a^2n(n+1)e^{-at}V_n^*(t)=0 \tag{2-5} \]

and form on the interval \((0,\infty)\) a complete orthogonal system with weight

\[ \frac{1}{1-e^{-at}}, \]

i.e.,

\[ \int_0^\infty \frac{V_n^*(t)V_m^*(t)}{1-e^{-at}}\,dt= \begin{cases} 0, & n\ne m,\\[4pt] \dfrac{1}{n(n+1)(2n+1)a}, & n=m, \end{cases} \tag{2-6} \]

with \(V_n^*(0)=V_n^*(\infty)=0\).

An explicit expression for the “e”-polynomials \(V_n^*(t)\) can be obtained by expanding the hypergeometric series in (2-3),

\[ V_n^*(t)=\frac{e^{-at}(1-e^{-at})}{n(n+1)} \sum_{k=0}^{n-1}(-1)^k \frac{(n+k+1)!}{(n-k-1)!k!(k+1)!}\,e^{-kat}, \tag{2-7} \]

For the first several values of “\(n\)” we shall have

\[ \left. \begin{aligned} V_1^*(t)&=e^{-at}(1-e^{-at}),\\ V_2^*(t)&=e^{-at}(1-e^{-at})(1-2e^{-at}),\\ V_3^*(t)&=e^{-at}(1-e^{-at})(1-5e^{-at}+5e^{-2at}),\\ V_4^*(t)&=e^{-at}(1-e^{-at})(1-9e^{-at}+21e^{-2at}-14e^{-3at}),\\ V_5^*(t)&=e^{-at}(1-e^{-at})(1-14e^{-at}+56e^{-2at}-84e^{-3at}+42e^{-4at}) \end{aligned} \right\} \tag{2-8} \]

Find the Laplace transform of the polynomial \(V_n^*(t)\). Introduce the notation: \(V_n(p)\doteq V_n^*(t)\) and \(P_n(p)\doteq P_n^*(t)\). According to (2-2) we can write

\[ 2(2n+1)V_n(p)=P_{n-1}(p)-P_{n+1}(p)\qquad (n=1,2,3,\ldots). \tag{2-9} \]

We have [4]

\[ P_{n-1}(p)=\frac{1}{p}\prod_{k=0}^{n-2}\frac{(k+1)a-p}{(k+1)a+p}, \]

\[ P_{n+1}(p)=\frac{1}{p}\prod_{k=0}^{n}\frac{(k+1)a-p}{(k+1)a+p}. \]

Formula (2-9), after transformations, gives

\[ V_n(p)=a\frac{\displaystyle\prod_{k=0}^{n-2}[(k+1)a-p]} {\displaystyle\prod_{k=0}^{n}[(k+1)a+p]} = \frac{a(a-p)(2a-p)\cdots[(n-1)a-p]} {(a+p)(2a+p)\cdots[(n+1)a+p]}, \]

where

\[ \prod_{k=0}^{n-2}[(k+1)a-p]=1 \]

for \(n=1\).

Let us consider the question of expanding an arbitrary function \(f(t)\) in a series in the integral polynomials \(V_n^*(t)\). Suppose that

\[ f(0)=f(\infty)=0. \tag{2-11} \]

By virtue of the completeness of the system \(V_n^*(t)\) \((n=1,2,\ldots)\), we can write

\[ f(t)=\sum_{n=1}^{\infty} C_n V_n^*(t). \tag{2-12} \]

Multiply both sides of (2-12) by

\[ \frac{V_m^*(t)}{1-e^{-at}} \]

and integrate over the interval \((0,\infty)\). Taking into account the properties (2-6), we shall have

\[ C_n=n(n+1)(2n+1)a\int_{0}^{\infty} \frac{f(t)V_n^*(t)}{1-e^{-at}}\,dt, \]

or, substituting expression (2-4),

\[ C_n=(2n+1)\int_{0}^{\infty} f(t)\frac{dP_n^*(t)}{dt}\,dt. \tag{2-13} \]

Finally, if the expanded representation (2-7) is taken into account, we obtain

\[ \frac{C_n}{n(n+1)(2n+1)a} = \sum_{k=0}^{n-1}\eta_{nk}M_k^* \quad (n=1,2,\ldots), \tag{2-14} \]

where \(M_k^*\) is the exponential moment of order \(k\) of the function \(f(t)\), and

\[ \eta_{nk}=(-1)^k \frac{(n+k+1)!}{(n-k-1)!\,k!\,(k+1)!\,n(n+1)}. \tag{2-15} \]

Or, in matrix-vector form,

\[ \overline{C}=\eta\overline{M}^{*}, \tag{2-16} \]

where \(\eta\) is a triangular matrix with components (2-15), with \(\eta_{nk}=0\) for \(k>n-1\), and \(\overline{C}\) and \(\overline{M}^{*}\) are the vectors of coefficients and moments, respectively, i.e.

\[ \overline{C}= \left\{ \frac{C_1}{6a},\, \frac{C_2}{30a},\,\ldots,\, \frac{C_n}{n(n+1)(2n+1)a},\,\ldots \right\}, \]

\[ \overline{M}^{*}=\{M_0^*,\,M_1^*,\,\ldots,\,M_{n-1}^*,\,\ldots\}. \]

In the Appendix the matrix \(\eta\) is given up to and including \(n=8\). Co-

The Fourier coefficients in the expansion (2-12) can be found by interpolation in the domain of operator images, analogously to how this was done for the coefficients of the expansion (1-2).

Let us note that if the function \(\varphi(t)\), to be determined from its operator image \(\varphi(p)\), does not vanish at the endpoints of the interval, i.e. \(\varphi(0)\ne 0\) and \(\varphi(\infty)\ne 0\), then it is sought in the form

\[ \varphi(t)=\varphi(0)e^{-at}+\varphi(\infty)(1-e^{-at})+\sum_{n=1}^{\infty} C_n V_n^*(t). \tag{2-17} \]

The expansion coefficients \(C_n\) \((n=1,2,\ldots)\) in this case are determined for the function

\[ f(t)=\varphi(t)-\varphi(0)e^{-at}-\varphi(\infty)(1-e^{-at}), \]

whose Laplace image has the form

\[ F(p)=\varphi(p)-\frac{\varphi(0)}{p+a}-\frac{a\varphi(\infty)}{p(p+a)}. \]

In other words, the exponential moments in this case must be determined by the formula

\[ M_{n-1}^* = F(na)=\varphi(na)-\frac{\varphi(0)}{(n+1)a} -\frac{\varphi(\infty)}{n(n+1)a}\quad (n=1,2,\ldots)^3 \tag{2-18} \]

and \(\varphi(0)\) and \(\varphi(\infty)\) can be found from the theorems on the initial and limiting values.

In practice, for purposes of inverting the Laplace transform it is more convenient to use the normalized integral Legendre “\(e\)” polynomials \(\overline{V}_n^*(t)\) [3]:

\[ \overline{V}_n^*(t)= \frac{2^n(n+1)\left[\Gamma\left(\frac{n+1}{2}\right)\right]^2}{(n-1)!}\,V_n^*(t). \tag{2-19} \]

For these polynomials the estimate holds

\[ \left|\overline{V}_n^*(t)\right| \leq 1 \quad 0\leq t\leq \infty. \]

If the expansion in the “\(e\)” polynomials \(V_n^*(t)\) is represented in the form

\[ f(t)\approx \sum_{n=1}^{m} A_n \overline{V}_n^*(t), \tag{2-20} \]

then the coefficients \(A_n\) will be related to the coefficients \(C_n\) by the following formula:

\[ A_n= \frac{(n-1)!}{2^n(n+1)\left[\Gamma\left(\frac{n+1}{2}\right)\right]^2}\,C_n. \tag{2-21} \]

Taking the Laplace transform of (2-20), and taking into account (2-10) and (2-19), we obtain

\[ F(p)\approx \sum_{n=1}^{m} A_n \frac{2^n(n+1)\left[\Gamma\left(\frac{n+1}{2}\right)\right]^2 a}{(n-1)!} \frac{\displaystyle\prod_{k=0}^{\,n-2}[(k+1)a-p]} {\displaystyle\prod_{k=0}^{\,n}[(k+1)a+p]}, \]

i.e. the operator image \(F(p)\) is approximated by a finite expansion, each term of which is a rational fraction.

3. Determination of the Parameter “\(a\)”

Up to now we have not touched at all on the question of the magnitude of the material parameter “\(a\).” It can be rigorously shown that infinite series of the form (1-2) or (2-20) will converge and represent the desired function for any finite value of the parameter “\(a\).” The inversion of the Laplace transform by means of exponential functions for \(a=1\) was considered in [1, 2]. Meanwhile, the value of the parameter “\(a\)” has a substantial influence on the accuracy of the approximation of the desired original by a segment of the series in the functions \(S_k(t)\) or the “\(\varepsilon\)”-polynomials \(V_n^*(t)\). Therefore the value \(a=1\) cannot, in general, provide the required accuracy in the approximation of the original by a finite number of terms of the series. For this, a practically unacceptable number of terms would be required. The series (1-2) and (2-12), asymptotically (i.e., for large values of \(t\)) and up to a constant factor, tend to zero as the exponentials

\[ e^{-\frac{at}{2}} \]

or \(e^{-at}\), respectively; consequently, the parameter “\(a\)” must be determined in such a way that the rate of decay of the original \(f(t)\) for large values of \(t\) is approximately equal to the rate of decay of the exponential \(e^{-\frac{at}{2}}\) or \(e^{-at}\).

If the Laplace transform is given in the form of a rational fraction, then the corresponding original is represented as a sum of exponentials with complex exponents in the general case. The rate of decay of such an original is determined by the term whose exponent has the smallest real part \((\alpha)\). This exponent can be isolated by the method described in [5]. The parameter “\(a\)” in this case is, obviously, equal to \(2\alpha\) or \(\alpha\), depending on which functions are used to represent the original. The essence of the method mentioned consists in the fact that the moments of the original, i.e., the quantities

\[ m_k=\int_0^\infty t^k f(t)\,dt \qquad (k=0,1,2,\ldots), \]

as \(k\) increases, are determined to an ever greater extent by the finite portion of the curve \(f(t)\), and the character of the variation of the moment sequence makes it possible to judge the character of the curve \(f(t)\) for large values of \(t\) [5]. This idea can be successfully applied to determine the parameter “\(a\)” in the case when the given image is not a rational fraction. Let the given image \(F(p)\) of some original \(f(t)\) be an analytic function in the entire right half-plane, including the imaginary axis; then all derivatives exist in neighborhoods of the origin. Suppose that they are all nonzero.

In this case

\[ \lim_{p\to 0} pF(p)=\lim_{t\to\infty} f(t)=0. \]

Taking into account that

\[ F^{(k)}(0)=(-1)^k\int_0^\infty t^k f(t)\,dt=(-1)^k m_k. \tag{3-1} \]

We can find the sequence of moments \(m_k\) \((k=1,2,\ldots,n)\) or the sequence of quantities

\[ M_k=\frac{1}{k!}\,m_k=\frac{(-1)^k}{k!}F^{(k)}(0). \]

If this sequence varies monotonically without changing sign, then this means that the original \(f(t)\), for large values of \(t\), tends to zero approximately according to the exponential law \(e^{-\alpha t}\). In this case, for sufficiently large values of \(k\)

\[ a \approx \frac{M_k}{M_{k+1}} . \tag{3-2} \]

If, however, the moment sequence \(M_k\) \((k=1,2,\ldots)\) changes sign, then the tendency of the original to zero is oscillatory in character, and the envelope of this oscillatory process varies approximately according to an exponential law. In this case [5]

\[ a \approx \frac{1}{2}\, \frac{\left|M_{k+2}\cdot M_{k-1}-M_{k+1}\cdot M_k\right|} {M_{k+1}^{2}-M_k\cdot M_{k+2}} . \tag{3-3} \]

Thus, \(a=2\alpha\) if the \(e\)-functions \(S_n(t)\) are used for the purposes of inversion, and \(a=\alpha\) if the \(e\)-polynomials \(V_n^*(t)\) are applied.

Sometimes the value of \(a\) determined by formula (3-2) does not stabilize as \(k\) increases. This means that, for large values of \(t\), the desired original differs appreciably from an exponential. In this case the parameter \(a\) can be found by the following device. In series (1-17) we restrict ourselves to only two terms, i.e., we set

\[ F(p) \approx B_1 S_1(p) + B_2 S_2(p). \tag{3-4} \]

For \(p=0\) we shall have

\[ F(0) \approx \frac{\pi}{a}(B_1-B_2). \]

Taking into account that

\[ B_1=\frac{4a}{\pi}F(a) \quad \text{and} \quad B_2=\frac{4a}{\pi}\{-2F(a)+4F(2a)\}, \]

we obtain

\[ \frac{1}{4}F(0)=3F(a)-4F(2a). \tag{3-5} \]

Solving equation (3-5) by any approximate method, we find a certain value of \(a\), which should then be adopted for the computation.

The idea underlying this device is that the behavior of the original for large values of \(t\) is determined by the behavior of its Laplace transform in neighborhoods of the origin; therefore, by determining the parameters of the segment of the expansion (3-4) from the condition of equality to the values of \(F(p)\) at \(p=0,a,2a\), we form, in the main, those values of the original that correspond to large values of \(t\).

Let us note that the value of \(a\) found should be reduced by a factor of two if the integral Legendre \(e\)-polynomials \(V_n^*(t)\) are used for the purposes of inversion. As an example, let us consider the inversion of a Laplace transform that is a transcendental function of the variable \(p\).

Let the operational image of some function \(\varphi(t)\) have the form

\[ \varphi(p)= \frac{\exp\left(\dfrac{b}{p+\beta}\right)} {p+\beta}. \tag{1-1} \]

Let us find an approximate expression for the original. We write (1-1) in the form

\[ (p+\beta)\varphi(p)=\exp\left(\frac{b}{p+\beta}\right). \]

Differentiating this expression with account of (1-1) gives the identity

\[ (p+b+\beta)\varphi(p)+(p+\beta)^2\varphi'(p)=0. \tag{1-II} \]

Differentiate it \(k\) times and put \(p=0\). Taking into account that

\[ \frac{d^k}{dp^k}(p+b+\beta)\varphi(p) =(p+b+\beta)\varphi^{(k)}(p)+k\varphi^{(k-1)}(p), \]

\[ \frac{d^k}{dp^k}(p+\beta)^2\varphi'(p) =(p+\beta)^2\varphi^{(k+1)}(p)+2k(p+\beta)\varphi^{(k)}(p)+k(k-1)\varphi^{(k-1)}(p), \]

we obtain the following recurrence relation for the derivatives at the origin:

\[ k^2\varphi^{(k-1)}(0)+[(2k+1)\beta+b]\varphi^{(k)}(0)+\beta^2\varphi^{(k+1)}(0)=0. \tag{1-III} \]

For definiteness put \(b=\beta=1\); then the relation takes the form

\[ k^2\varphi^{(k-1)}(0)+2(k+1)\varphi^{(k)}(0)+\varphi^{(k+1)}(0)=0 \quad (k=0,1,2,\ldots). \]

Successively we easily find \(\varphi'(0), \varphi''(0), \ldots\), and then also the quantities

\[ M_k=\frac{(-1)^k\varphi^{(k)}(0)}{k!}. \]

The computations give:

\[ \begin{gathered} M_1=2\varphi(0);\quad M_2=3.5\varphi(0);\quad M_3=8.5\varphi(0);\quad M_4=8.7\varphi(0);\quad M_5=12.8\varphi(0);\\ M_6=18.5\varphi(0);\quad M_7=26\varphi(0);\quad M_8=35.9\varphi(0);\quad M_9=48.5\varphi(0);\quad M_{10}=64.8\varphi(0). \end{gathered} \]

The moment sequence changes monotonically without a change of sign; consequently, we can use formula (3-2). For various values of \(k\) the formula gives

\[ \frac{M_4}{M_5}=0.68;\quad \frac{M_5}{M_9}=0.69;\quad \frac{M_6}{M_7}=0.71;\quad \frac{M_7}{M_8}=0.725. \]

\[ \frac{M_8}{M_9}=0.74;\quad \frac{M_9}{M_{10}}=0.75. \]

Thus, as \(k\) increases, the value of \(\alpha\) does not stabilize. Let us find “\(a\)” from equation (3-5). Since

\[ \lim_{p\to\infty} p\varphi(p)=1, \]

we shall have

\[ F(p)=\varphi(p)-\frac{1}{p+\frac{a}{2}} =\frac{1}{p+1}e^{\frac{1}{p+1}}-\frac{2}{2p+a}, \]

whence we easily obtain

\[ F(0)=e-\frac{2}{a}, \]

\[ F(a)=\frac{1}{a+1}e^{\frac{1}{a+1}}-\frac{2}{3a}, \]

\[ F(2a)=\frac{1}{2a+1}e^{\frac{1}{2a+1}}-\frac{2}{5a}. \]

Equation (3-5) takes the form

\[ \frac{1}{4}e-\frac{1}{10a} = \frac{3}{a+1}e^{\frac{1}{a+1}} -\frac{4}{2a+1}e^{\frac{1}{2a+1}}. \]

The approximate value of the positive root of this equation, i.e., of the parameter “\(a\),” is equal to \(1{,}09\) (the calculation was carried out on a slide rule). The calculation of the first 6 exponential moments is given in Table 1. The values of the coefficients \(B_n\) for the expansion in “\(e\)” functions of Chebyshev of the third kind \(S_n(t)\) are also given there. If we restrict ourselves to 4 terms, then the sought original is written in the form

\[ \varphi(t) \approx e^{-0{,}545t} + 0{,}22266 S_1(t) - 0{,}0917 S_2(t) + 0{,}00565 S_3(t) + 0{,}0054 S_4(t), \]

where the maximum absolute error will be no greater than \(0{,}0015\).

Table 1

Table 2

Table 2 gives the results of calculating the coefficients for representing the original in the form of a segment of a series in normalized integral “\(e\)” polynomials \(\overline{V}_n^*(t)\). The value of the parameter “\(a\)” in this case is equal to \(0{,}545\). The approximate representation has the form

\[ \varphi(t) \approx e^{-0{,}545t} + 0{,}26156\,\overline{V}_1^*(t) - 0{,}00325\,\overline{V}_2^*(t) - 0{,}00974\,\overline{V}_3^*(t) - \]

\[ -0{,}00708\,\overline{V}_4^*(t) - 0{,}0044\,\overline{V}_5^*(t). \]

In conclusion, we note that the exact value of the original is equal to

\[ \varphi(t)=e^{-t}I_0(2\sqrt{t}). \]

where \(I_0(2\sqrt{t})\) is the modified cylindrical function of zero order. We also note that the computation of the exponential moments \(M^*_{n-1}\) must be carried out with a very high degree of accuracy because of the rapid growth of the elements of the matrices \(\beta\) and \(\eta\), whereas for determining the parameter “\(a\)” the accuracy of an ordinary slide rule is sufficient.

Appendix II

1. Matrix \(\beta\)

\[ \beta_{kn} = \frac{(-1)^{n+k-1} 2^n (n+k)!} {(2n+1)(k-n-1)!\,n!\,(2n-1)!!} \]

2. Matrix \(\eta\)

\[ \eta_{nk} = (-1)^k \frac{(n+k+1)!} {(n-k-1)!\,k!\,(k+1)!\,n(n+1)} \]

REFERENCES

1. G. Doetsch, Guide to the Practical Application of the Laplace Transform. Nauka, Moscow, 1965.
2. K. Lanczos, Practical Methods of Applied Analysis. Fizmatgiz, Moscow, 1961.
3. V. M. Osipov, “Integral Legendre polynomials and the approximation of functions.” In the present volume.
4. V. M. Osipov, “Exponential polynomials and the expansion of certain typical signals.” Izv. TPI, vol. 180, Tomsk, 1969.
5. V. M. Osipov, “Determination of the zeros and poles of a minimum-phase transfer function by the method of moment sequences.” Izv. TPI, vol. 192, Tomsk, 1969.