InverseLaplace

ft::Talbot = Talbot(func::Function, Nterms::Integer=32)

Return ft, which estimates the inverse Laplace transform of func with the fixed Talbot algorithm. ft(t) evaluates the transform at t. You may want to tune Nterms together with setprecision(BigFloat, x).

Example

Compute the inverse transform of the transform of cos at argument pi/2.

julia> ft = Talbot(s -> s / (s^2 + 1), 80);

julia> ft(pi / 2)
-3.5366510684573195e-5

Note that given Float64 input, the precision of the returned value may not be satisfactory.

julia> Float64(ft(big(pi) / 2))
2.114425886215604e-49

ILT(function, Nterms=32)

This is an alias for the default Talbot() method.

ft::GWR = GWR(func::Function, Nterms::Integer=16)

Return ft, which estimates the inverse Laplace transform of func with the GWR algorithm. ft(t) evaluates the transform at t.

Example

Compute the inverse transform of the transform of cos at argument pi / 2.

julia> ft = GWR(s -> s / (s^2 + 1), 16);

julia> ft(pi / 2)
-0.001

w::Weeks{datatype} = Weeks(func::Function, Nterms::Integer=64, sigma=1.0, b=1.0; datatype=Float64)

Return w, which estimates the inverse Laplace transform of func with the Weeks algorithm. w(t) evaluates the transform at t. The accuracy depends on the choice of sigma and b, with the optimal choices depending on t. datatype should agree with the DataType returned by func. For convenience, datatype=Complex is equivalent to datatype=Complex{Float64}

The call to Weeks that creates w is expensive relative to evaluation via w(t).

Example

Compute the inverse transform of the transform of cos at argument pi/2.

julia> ft = Weeks(s -> s/(s^2+1), 80);

julia> ft(pi/2)
0.0

w::WeeksErr{datatype} = WeeksErr(func::Function, Nterms::Integer=64, sigma=1.0, b=1.0; datatype=Float64)

Return w, which estimates the inverse Laplace transform of func via the Weeks algorithm. w(t) returns a tuple containing the inverse transform at t and an error estimate. The accuracy of the inversion depends on the choice of sigma and b. See the documentation for Weeks for a description of the parameter datatype.

Example

Compute the inverse transform of the transform of cos, and an error estimate, at argument pi/2 using 80 terms.

julia> ft = WeeksErr(s -> s/(s^2+1), 80);

julia> ft(pi/2)
(0.0,3.0872097665938698e-15)

This estimate is more accurate than cos(pi/2).

julia> ft(pi/2)[1] - cos(pi/2)
-6.123233995736766e-17

julia> ft(pi/2)[1] - 0.0         # exact value
0.0

julia> ft(pi/2)[1] - cospi(1/2)  # cospi is more accurate
0.0

gaverstehfest(func::Function, t::AbstractFloat, v = stehfest_coeffs(N))

Evaluate the inverse Laplace transform of func at the point t using the Gaver-Stehfest algorithm. v are the stehfest coefficients computed with stehfest_coeffs that only depend on number of terms. N (which must be even) defaults to 36 should depend on the precision used and desired accuracy. The precision of the stehfest coefficients is used in the computation.

In double precision (Float64), N should be <= 18 to provide best accuracy. Usually only moderate accuracy can be achieved in double precision ~rtol ≈ 1e-4. For higher accuracy, BigFloats should be used with N = 36.

Increasing precision should be accompanied by an increase in the number of coefficients used. Increasing precision without increasing number of coefficients will not yield better accuracy. The inverse is generally true as well.

This method is not robust to oscillating F(t) and must be smooth.

Examples

julia> F(s) = 1 / (s + 1) # where analytical inverse is f(t) = exp(-t)

julia> InverseLaplace.gaverstehfest(F, 2.0) # computes with default 36 coefficients
0.1353352832366128315426471959891170863784520272858265151734040491181311503059561

julia> InverseLaplace.gaverstehfest(F, 2.0, v = stehfest_coeffs(20)) # computes with custom number
0.1353353114073885136885007645878189977740624364169818512599342310325432753817199

# to calculate in double precision convert stehfest coefficients to Float64
julia> InverseLaplace.gaverstehfest(F, 2.0, v = convert(Vector{Float64}, stehfest_coeffs(18)))
0.13533650985980258

gaverstehfest(func::Function, t::AbstractArray, v = stehfest_coeffs(N))

Evaluate the inverse Laplace transform of func over an array of t using the Gaver-Stehfest algorithm. Computes coefficients once and calculates f(t) across available threads.

Example

julia> gaverstehfest(s -> 1/(s + 1), 2.0:4.0)
0.1353355048178631463198819259249043857320738467297190227379428405365703748720082
0.04978728177016550841951683309410878070827215175940039571536203126773133604403101
0.01831383641619355549133471411401755204406266755309342902744499924574072011058623

setNterms(ailt::AbstractILt, Nterms::Integer)

set the number of terms used in the inverse Laplace tranform ailt. If ailt stores internal data, it will be recomputed, so that subsequent calls ailt(t) reflect the new value of Nterms.

setparameters(w::AbstractWeeks, sigma, b, Nterms)

Set the parameters for the inverse Laplace transform object w and recompute the internal data. Subsequent calls w(t) will use these parameters. If Nterms or both Nterms and b are omitted, then their current values are retained.

ILtPair(ilt::AbstractILt, ft::Function)

return an object of type ILtPair that associates ilt the inverse Laplace transform of a function with it "exact" numerical inverse ft. Calling abserr(p, t) returns the absolute error between the inverse transform and the exact value.

Example

This example compares the inversion using the Weeks algorithm of the Laplace transform of cos(t) to its exact value at t=1.0.

julia> p = ILtPair(Weeks(s -> s/(1+s^2)), cos);
julia> abserr(p, 1.0)

0.0

abserr(p::ILtPair, t)

Compute the absolute error between the estimated inverse Laplace transform and "exact" numerical solution contained in p at the point t. See ILtPair.

iltpair_power(n)

Return a TransformPair for the power function x^n. This can be used to construct an ILTPair.

julia> p  = Talbot(iltpair_power(3));

julia> Float64(abserr(p, 1)) # test Talbot method for Laplace transform of x^3, evaluated at 1
2.0820366247539812e-26

TransformPair{T,V}

A pair of functions for analyzing an inverse Laplace transform method. Field ft is the real-space function. Field fs is the Laplace-space function.

ilt(func::Function, t::AbstractFloat, M::Integer=32)

ilt is an alias for the default inverse Laplace transform method talbot.

talbot(func::Function, t::AbstractFloat, M::Integer=talbot_default_num_terms)

Evaluate the inverse Laplace transform of func at the point t. Use M terms in the algorithm. For typeof(t) is Float64, the default for M is 32. For BigFloat the default is 64.

If BigFloat precision is larger than default, try increasing M.

Example

julia> InverseLaplace.talbot(s -> 1 / s^3, 3)
4.50000000000153

gwr(func::Function, t::AbstractFloat, M::Integer=gwr_default_num_terms)

Evaluate the inverse Laplace transform of func at the point t. Use M terms in the algorithm. For typeof(t) is Float64, the default for M is 16. For BigFloat the default is 64.

If BigFloat precision is larger than default, try increasing M.

Example

julia> InverseLaplace.gwr( s -> 1/s^3, 3.0)
4.499985907607361

talbotarr(func, ta::AbstractArray, M)

Compute the inverse Laplace transform for each element in ta. Each evaluation of func(s) is used for all elements of ta. This may be faster than a broadcast application of talbot (i.e. talbot.(...) , but is in general, less accurate. talbotarr uses the "fixed" Talbot method.