Inverse Laplacian
Numerical inversion of Laplace transforms.
invlap(Fs, t1, t2, nnt, a = 6, ns = 20, nd = 19)
Fs |
function representing the function to be inverse-transformed. |
t1, t2 |
end points of the interval. |
nnt |
number of grid points between t1 and t2. |
a |
shift parameter; it is recommended to preserve value 6. |
ns, nd |
further parameters, increasing them leads to lower error. |
The transform Fs may be any reasonable function of a variable s^a, where a
is a real exponent. Thus, the function invlap can solve fractional
problems and invert functions Fs containing (ir)rational or transcendental
expressions.
Returns a list with components x the x-coordinates and y
the y-coordinates representing the original function in the interval
[t1,t2].
Based on a presentation in the first reference. The function invlap
on MatlabCentral (by ) served as guide. The Talbot procedure from the
second reference could be an interesting alternative.
J. Valsa and L. Brancik (1998). Approximate Formulae for Numerical Inversion of Laplace Transforms. Intern. Journal of Numerical Modelling: Electronic Networks, Devices and Fields, Vol. 11, (1998), pp. 153-166.
L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer (2006). Talbot quadratures and rational approximations. BIT. Numerical Mathematics, 46(3):653–670.
Fs <- function(s) 1/(s^2 + 1) # sine function Li <- invlap(Fs, 0, 2*pi, 100) ## Not run: plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid() Fs <- function(s) tanh(s)/s # step function L1 <- invlap(Fs, 0.01, 20, 1000) plot(L1[[1]], L1[[2]], type = "l", col = "blue") L2 <- invlap(Fs, 0.01, 20, 2000, 6, 280, 59) lines(L2[[1]], L2[[2]], col="darkred"); grid() Fs <- function(s) 1/(sqrt(s)*s) L1 <- invlap(Fs, 0.01, 5, 200, 6, 40, 20) plot(L1[[1]], L1[[2]], type = "l", col = "blue"); grid() Fs <- function(s) 1/(s^2 - 1) # hyperbolic sine function Li <- invlap(Fs, 0, 2*pi, 100) plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid() Fs <- function(s) 1/s/(s + 1) # exponential approach Li <- invlap(Fs, 0, 2*pi, 100) plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid() gamma <- 0.577215664901532 # Euler-Mascheroni constant Fs <- function(s) -1/s * (log(s)+gamma) # natural logarithm Li <- invlap(Fs, 0, 2*pi, 100) plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid() Fs <- function(s) (20.5+3.7343*s^1.15)/(21.5+3.7343*s^1.15+0.8*s^2.2+0.5*s^0.9)/s L1 <- invlap(Fs, 0.01, 5, 200, 6, 40, 20) plot(L1[[1]], L1[[2]], type = "l", col = "blue") grid() ## End(Not run)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.