SINLP
Computes the inverse Laplace transform of a smooth complex function.
Required Arguments
F — User-supplied FUNCTION to which the inverse Laplace transform will be
computed. The form is F(Z), where
Z — Complex argument. (Input)
F — The complex function value. (Output)
F must be declared EXTERNAL in the calling program. F must also be declared COMPLEX.
T — Vector of length N containing points at which the inverse Laplace transform is desired. (Input)
T(I) must be greater than zero for all I.
FINV — Vector of length N whose I-th component contains the approximate value of the inverse Laplace transform at the point T(I). (Output)
Optional Arguments
N — The number of points at which the inverse Laplace transform is desired. (Input)
Default: N = size (T,1).
SIGMA0 — An estimate for the maximum of the real parts of the singularities of F. (Input)
If unknown, set SIGMA0 = 0.0.
Default: SIGMA0 = 0.e0.
EPSTOL — The required absolute uniform pseudo accuracy for the coefficients and inverse Laplace transform values. (Input)
Default: EPSTOL = 1.1920929e-5 for single precision and 2.22d-10 for double precision.
ERRVEC — Vector of length eight containing diagnostic information. (Output)
All components depend on the intermediately generated Laguerre coefficients. See
Comments.
FORTRAN 90 Interface
Generic: CALL SINLP (F, T, FINV [, …])
Specific: The specific interface names are S_SINLP and D_SINLP.
FORTRAN 77 Interface
Single: CALL SINLP (F, N, T, SIGMA0, EPSTOL, ERRVEC, FINV)
Double: The double precision name is DSINLP.
Description
The routine SINLP computes the inverse Laplace transform of a complex-valued function. Recall that if f is a function that vanishes on the negative real axis, then we can define the Laplace transform of f by
It is assumed that for some value of s the integrand is absolutely integrable.
The computation of the inverse Laplace transform is based on a modification of Weeks’ method (see W.T. Weeks (1966)) due to B.S. Garbow et. al. (1988). This method is suitable when f has continuous derivatives of all orders on [0, ∞). In this situation, this routine should be used in place of the IMSL routine INLAP. It is especially efficient when multiple function values are desired. In particular, given a complex-valued function F(s) = L[f ](s), we can expand f in a Laguerre series whose coefficients are determined by F. This is fully described in B.S. Garbow et. al. (1988) and Lyness and Giunta (1986).
The algorithm attempts to return approximations g(t) to f(t) satisfying
where ɛ := EPSTOL and σ := SIGMA > SIGMA0. The expression on the left is called the pseudo error. An estimate of the pseudo error is available in ERRVEC(1).
The first step in the method is to transform F to ɸ where
Then, if f is smooth, it is known that ɸ is analytic in the unit disc of the complex plane and hence has a Taylor series expansion
which converges for all z whose absolute value is less than the radius of convergence Rc. This number is estimated in ERRVEC(6). In ERRVEC(5), we estimate the smallest number K which satisfies
for all R < Rc.
The coefficients of the Taylor series for ɸ can be used to expand f in a Laguerre series
Comments
1. Workspace may be explicitly provided, if desired, by use of S2NLP/DS2NLP. The reference is:
CALL S2NLP (F, N, T, SIGMA0, EPSTOL, ERRVEC, FINV, SIGMA, BVALUE, MTOP, WK, IFLOVC)
The additional arguments are as follows:
SIGMA — The first parameter of the Laguerre expansion. If SIGMA is not greater than SIGMA0, it is reset to SIGMA0 + 0.7. (Input)
BVALUE — The second parameter of the Laguerre expansion. If BVALUE is less than 2.0 * (SIGMA ‑ SIGMA0), it is reset to 2.5 * (SIGMA ‑ SIGMA0). (Input)
MTOP — An upper limit on the number of coefficients to be computed in the Laguerre expansion. MTOP must be a multiple of four. Note that the maximum number of Laplace transform evaluations is MTOP/2 + 2. (Default: 1024.) (Input)
WK — Real work vector of length 9 * MTOP/4.
IFLOVC — Integer vector of length N, the I-th component of which contains the overflow/underflow indicator for the computed value of FINV(I). (Output) See Comment 3.
2. Informational errors
Type | Code | Description |
|---|
1 | 1 | Normal termination, but with estimated error bounds slightly larger than EPSTOL. Note, however, that the actual errors on the final results may be smaller than EPSTOL as bounds independent of T are pessimistic. |
3 | 2 | Normal calculation, terminated early at the roundoff error level estimate because this estimate exceeds the required accuracy (usually due to overly optimistic expectation by the user about attainable accuracy). |
4 | 3 | The decay rate of the coefficients is too small. It may improve results to use S2NLP and increase MTOP. |
4 | 4 | The decay rate of the coefficients is too small. In addition, the roundoff error level is such that required accuracy cannot be reached. |
4 | 5 | No error bounds are returned as the behavior of the coefficients does not enable reasonable prediction. Results are probably wrong. Check the value of SIGMA0. In this case, each of ERRVEC(J), J = 1, …, 5, is set to -1.0. |
3. The following are descriptions of the vectors ERRVEC and IFLOVC.
ERRVEC — Real vector of length eight.
ERRVEC(1) = Overall estimate of the pseudo error, ERRVEC(2) + ERRVEC(3) + ERRVEC(4). Pseudo error = absolute error / exp(sigma * tvalue).
ERRVEC(2) = Estimate of the pseudo discretization error.
ERRVEC(3) = Estimate of the pseudo truncation error.
ERRVEC(4) = Estimate of the pseudo condition error on the basis of minimal noise levels in the function values.
ERRVEC(5) = K, the coefficient of the decay function for ACOEF, the coefficients of the Laguerre expansion.
ERRVEC(6) = R, the base of the decay function for ACOEF. Here
abs(ACOEF (J + 1)).LE.K/R**J for J.GE.MACT/2, where MACT is the number of Laguerre coefficients actually computed.
ERRVEC(7) = ALPHA, the logarithm of the largest ACOEF.
ERRVEC(8) = BETA, the logarithm of the smallest nonzero ACOEF.
IFLOVC — Integer vector of length N containing the overflow/underflow indicators for FINV. For each I, the value of IFLOVC(I) signifies the following.
0 = Normal termination.
1 = The value of the inverse Laplace transform is found to be too large to be representable; FINV(I) is set to AMACH(6).
–1 = The value of the inverse Laplace transform is found to be too small to be representable; FINV(I) is set to 0.0.
2 = The value of the inverse Laplace transform is estimated to be too large, even before the series expansion, to be representable; FINV(I) is set to AMACH(6).
–2 = The value of the inverse Laplace transform is estimated to be too small, even before the series expansion, to be representable; FINV(I) is set to 0.0.
Example
We invert the Laplace transform of the simple function (s - 1)-2 and print the computed answer, the true solution, and the difference at five different points. The correct inverse transform is xex.
USE SINLP_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER I, NOUT
REAL DIF(5), ERRVEC(8), EXP, FINV(5), FLOAT, RELERR, &
SIGMA0, T(5), TRUE(5), EPSTOL
COMPLEX F
INTRINSIC EXP, FLOAT
EXTERNAL F
! Get output unit number
CALL UMACH (2, NOUT)
!
DO 10 I=1, 5
T(I) = FLOAT(I) - 0.5
10 CONTINUE
SIGMA0 = 1.0E0
RELERR = 5.0E-4
EPSTOL = 1.0E-4
CALL SINLP (F, T, FINV, SIGMA0=SIGMA0, EPSTOL=RELERR)
! Evaluate the true solution and the
! difference
DO 20 I=1, 5
TRUE(I) = T(I)*EXP(T(I))
DIF(I) = TRUE(I) - FINV(I)
20 CONTINUE
!
WRITE (NOUT,99999) (T(I),FINV(I),TRUE(I),DIF(I),I=1,5)
99999 FORMAT (7X, 'T', 8X, 'FINV', 9X, 'TRUE', 9X, 'DIFF', /, &
5(1X,E9.1,3X,1PE10.3,3X,1PE10.3,3X,1PE10.3,/))
END
!
COMPLEX FUNCTION F (S)
COMPLEX S
!
F = 1./(S-1.)**2
RETURN
END
Output
T FINV TRUE DIFF
0.5E+00 8.244E-01 8.244E-01 -2.086E-06
1.5E+00 6.723E+00 6.723E+00 -8.583E-06
2.5E+00 3.046E+01 3.046E+01 0.000E+00
3.5E+00 1.159E+02 1.159E+02 2.289E-05
4.5E+00 4.051E+02 4.051E+02 -2.136E-04
