.p>.CMCH6.DOC!INVERSE_LAPLACE;inverse_laplace
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 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 e 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 e = pseudo_accuracy and s = 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 f where
Then, if f is smooth, it is known that f 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 f can be used to expand f in a Laguerre series
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>
main()
{
f_complex f(f_complex);
int n = 5;
float t[5];
float true_inverse[5];
float relative_diff[5];
int i;
float *inverse;
/* 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>
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |