Computes the inverse Laplace transform of a complex function.
Synopsis
#include <imsl.h>
float *imsl_f_inverse_laplace (f_complex fcn(), float sigma0, int n, float t[], …, 0)
The type double procedure is imsl_d_inverse_laplace.
Required Arguments
f_complex fcn(f_complex z)
(Input)
User-supplied function for which the inverse Laplace transform will
be computed.
float sigma0
(Input)
An estimate for the maximum of the real parts of the singularities of
fcn. If unknown,
set sigma0 =
0.0.
int n
(Input)
The number of points at which the inverse Laplace transform is
desired.
float t[]
(Input)
Array of size n containing the points at
which the inverse Laplace transform is desired.
Return Value
A pointer to the array of length n whose i-th component contains the approximate value of the inverse Laplace transform at the point t[i]. To release this space, use imsl_free. If no solution was computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
float *imsl_f_inverse_laplace (f_complex fcn(), float sigma0, int n, float t[],
IMSL_RETURN_USER, float x[],
IMSL_PSEUDO_ACCURACY, float pseudo_accuracy,
IMSL_FIRST_LAGUERRE_PARAMETER, float sigma,
IMSL_SECOND_LAGUERRE_PARAMETER, float bvalue,
IMSL_MAXIMUM_COEFFICIENTS, int mtop,
IMSL_ERROR_EST, float *error_est,
IMSL_DISCRETIZATION_ERROR_EST, float *disc_error_est,
IMSL_TRUNCATION_ERROR_EST, float *trunc_error_est,
IMSL_CONDITION_ERROR_EST, float *cond_error_est,
IMSL_DECAY_FUNCTION_COEFFICIENT, float *k,
IMSL_DECAY_FUNCTION_BASE, float *r,
IMSL_LOG_LARGEST_COEFFICIENTS, float *log_largest_coefs,
IMSL_LOG_SMALLEST_COEFFICIENTS, float *log_smallest_coefs,
IMSL_UNDER_OVERFLOW_INDICATORS, Imsl_laplace_flow **indicators,
IMSL_FCN_W_DATA, f_complex fcn(), void *data,
0)
Optional Arguments
IMSL_RETURN_USER, float x[]
(Output)
A user-allocated array of length n containing the approximate
value of the inverse Laplace transform.
IMSL_PSEUDO_ACCURACY, float pseudo_accuracy
(Input)
The required absolute uniform pseudo accuracy for the coefficients
and inverse Laplace transform values.
Default: pseudo_accuracy =
, where ε is machine epsilon
IMSL_FIRST_LAGUERRE_PARAMETER, float sigma
(Input)
The first parameter of the Laguerre expansion. If sigma is not greater
than sigma0, it
is reset to sigma0 + 0.7.
Default:
sigma = sigma0 + 0.7
IMSL_SECOND_LAGUERRE_PARAMETER, float bvalue
(Input)
The second parameter of the Laguerre expansion. If bvalue is less than
2.0*(sigma − sigma0), it is reset
to 2.5*(sigma − sigma0).
Default:
bvalue =
2.5*(sigma − sigma0)
IMSL_MAXIMUM_COEFFICIENTS, int mtop
(Input)
An upper limit on the number of coefficients to be computed in the
Laguerre expansion. Argument mtop must be a
multiple of four.
Default: mtop = 1024
IMSL_ERROR_EST, float *error_est
(Output)
Overall estimate of the pseudo error, disc_error_est + trunc_error_est + cond_error_est. See
the Description section for
details.
IMSL_DISCRETIZATION_ERROR_EST, float *disc_error_est
(Output)
Estimate of the pseudo discretization error.
IMSL_TRUNCATION_ERROR_EST, float *trunc_error_est
(Output)
Estimate of the pseudo truncation error.
IMSL_CONDITION_ERROR_EST, float *cond_error_est
(Output)
Estimate of the pseudo condition error on the basis of minimal noise
levels in the function values.
IMSL_DECAY_FUNCTION_COEFFICIENT, float *k
(Output)
The coefficient of the decay function. See the Description section for
details.
IMSL_DECAY_FUNCTION_BASE, float *r
(Output)
The base of the decay function. See the Description section for
details.
IMSL_LOG_LARGEST_COEFFICIENTS, float *log_largest_coefs
(Output)
The logarithm of the largest coefficient in the
decay function. See the Description section for
details.
IMSL_LOG_SMALLEST_COEFFICIENTS,
float *log_smallest_coefs
(Output)
The logarithm of the smallest nonzero coefficient in the
decay function. See the Description section for
details.
IMSL_UNDER_OVERFLOW_INDICATORS,
Imsl_laplace_flow **indicators
(Output)
The address of a pointer initialized by imsl_f_inverse_laplace
to point to an array of length n containing the
overflow/underflow indicators for the computed approximate inverse Laplace
transform. For the ith point at which the transform is computed,
indicators[i] signifies the following:
|
indicators [i] |
Meaning |
|
IMSL_NORMAL_TERMINATION |
Normal termination. |
|
IMSL_TOO_LARGE |
The value of the inverse Laplace transform is too large to be representable. This component of the result is set to NaN. |
|
IMSL_TOO_SMALL |
The value of the inverse Laplace transform is found to be too small to be representable. This component of the result is set to 0.0. |
|
IMSL_TOO_LARGE_BEFORE_EXPANSION |
The value of the inverse Laplace transform is estimated to be too large, even before the series expansion, to be representable. This component of the result is set to NaN. |
|
IMSL_TOO_SMALL_BEFORE_EXPANSON |
The value of the inverse Laplace transform is estimated to be too small, even before the series expansion, to be representable. This component of the result is set to 0.0. |
IMSL_FCN_W_DATA, f_complex fcn(f_complex z, void *data)
,void *data, (Input)
User
supplied function for which the inverse Laplace transform will be computed,
which also accepts a pointer to data that is supplied by the user. data is a pointer to
the data to be passed to the user-supplied function. See the Introduction, Passing Data to
User-Supplied Functions at the beginning of this manual for more
details.
Description
The function imsl_f_inverse_laplace 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 the Laplace transform of f is defined 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 Weeks (1966)) due to Garbow et al. (1988). This method is suitable when f has continuous derivatives of all orders on [0, ∞). In particular, given a complex-valued function F(s) = L[f] (s), f can be expanded in a Laguerre series whose coefficients are determined by F. This is fully described in Garbow et al. (1988) and Lyness and Giunta (1986).
The algorithm attempts to return approximations g(t) to f(t) satisfying
where ε = pseudo_accuracy and σ = sigma > sigma0. The expression on the left is called the pseudo error. An estimate of the pseudo error in available in error_est.
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 r, obtained through the optional argument IMSL_DECAY_FUNCTION_BASE. Using optional argument IMSL_DECAY_FUNCTION_COEFFICIENT, the smallest number K is estimated which satisfies
for all R < Rc.
The coefficients of the Taylor series for φ can be used to expand f in a Laguerre series
On some platforms, imsl_f_inverse_laplace can evaluate the user-supplied function fcn in parallel. This is done only if the function imsl_omp_options is called to flag user-defined functions as thread-safe. A function is thread-safe if there are no dependencies between calls. Such dependencies are usually the result of writing to global or static variables.
Examples
Example 1
This example computes the inverse Laplace transform of the function (s − 1)-2, and prints the computed approximation, true transform value, and difference at five points. The correct inverse transform is xex. From Abramowitz and Stegun (1964).
#include <imsl.h>
#include <math.h>
int main()
{
f_complex f(f_complex);
int n = 5;
float t[5];
float true_inverse[5];
float relative_diff[5];
int i;
float *inverse;
imsl_omp_options(IMSL_SET_FUNCTIONS_THREAD_SAFE, 1, 0);
/* Initialize t and compute inverse */
for (i=0; i<n; i++)
t[i] = (float)i + 0.5;
inverse = imsl_f_inverse_laplace(f, 1.5, n, t, 0);
/* Compute true inverse, relative difference */
for (i=0; i<n; i++) {
true_inverse[i] = t[i]*exp(t[i]);
relative_diff[i] = fabs(inverse[i] - true_inverse[i])/
true_inverse[i];
}
printf("\t t\t\t f_inv\t\t true\t\t diff\n");
for (i=0; i<n; i++)
printf ("\t%5.1f\t\t%7.3f\t\t%7.3f\t\t%6.1e\n", t[i],
inverse[i], true_inverse[i], relative_diff[i]);
}
f_complex f(f_complex s)
{
/* Return 1/(s-1)**2 */
f_complex one = {1.0, 0.0};
return (imsl_c_div(one,
imsl_c_mul(imsl_c_sub(s, one), imsl_c_sub(s, one))));
}
Output
t f_inv true diff
0.5 0.824 0.824 1.5e-05
1.5 6.722 6.723 1.0e-05
2.5 30.456 30.456 5.6e-07
3.5 115.906 115.904 1.8e-05
4.5 405.054 405.077 5.8e-05
Example 2
This example computes the inverse Laplace transform of the function e-1/s/s, and prints the computed approximation, true transform value, and difference at five points. Additionally, the inverse is returned in user-supplied space, and a required accuracy for the inverse transform values is specified. The correct inverse transform is
From Abramowitz and Stegun (1964).
#include <imsl.h>
#include <math.h>
int main()
{
f_complex f(f_complex);
int n = 5;
int i;
float t[5];
float true_inverse[5];
float relative_diff[5];
float inverse[5];
Imsl_laplace_flow *indicators;
/* Initialize t and compute inverse */
for (i=0; i<n; i++) t[i] = (float)i + 0.5;
imsl_f_inverse_laplace(f, 0.0, n, t,
IMSL_PSEUDO_ACCURACY, 1.0e-6,
IMSL_UNDER_OVERFLOW_INDICATORS, &indicators,
IMSL_RETURN_USER, inverse,
0);
/* Compute true inverse, relative
difference */
for (i=0; i<n; i++) {
true_inverse[i] = imsl_f_bessel_J0(2.0*sqrt(t[i]));
relative_diff[i] = fabs((inverse[i] - true_inverse[i])/
true_inverse[i]);
}
/* Print results, noting if any results
overflowed or underflowed */
printf("\t T\t\t f_inv\t\t true\t\t diff\n");
for (i=0; i<n; i++)
if (indicators[i] == IMSL_NORMAL_TERMINATION)
printf ("\t%5.1f\t\t%7.3f\t\t%7.3f\t\t%6.1e\n",
t[i],
inverse[i], true_inverse[i],
relative_diff[i]);
else
printf("Overflow or underflow noted.\n");
}
f_complex f(f_complex s)
{
/* Return (1/s)(exp(-1/s) */
f_complex one = {1.0, 0.0};
f_complex s_inverse;
s_inverse = imsl_c_div(one, s);
return (imsl_c_mul(s_inverse, imsl_c_exp(imsl_c_neg(s_inverse))));
}
Output
T f_inv true diff
0.5 0.559 0.559 2.1e-07
1.5 -0.023 -0.023 8.5e-06
2.5 -0.310 -0.310 9.6e-08
3.5 -0.401 -0.401 7.4e-08
