inverse_laplace

CNLMath : Transforms : inverse_laplace

inverse_laplace

Synopsis

Required Arguments

Return Value

Synopsis with Optional Arguments

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,

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 Descriptionsection 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_EXPANSION	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 Passing Data to User-Supplied Functions in the introduction to 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 <stdio.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 <stdio.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

4.5 -0.370 -0.370 6.4e-07

Fatal Errors

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm. User flag = "#".