Advanced Search
All Posts :
Here is a simple procedure doing numerical inverse Laplace transform using Post's inversion formula, 
ILT:=proc(f,s) local k,dig;
if nargs>2 and type(args[-1],'posint') then dig:=args[-1] 
else dig:=Digits fi;
proc(T) local t,d,a;
d:=dig;
a:=Limit();
t:=evalf(T,dig);
while op(0,a)='Limit' do
a:=evalf(Limit((-1)^k/k!*(k/t)^(k+1)*eval(diff(f,s$k),s=k/t),k=infinity),d); 
d:=2*d od;
evalf(a,dig) end end:
The original version of the procedure looked much more simple, 
ILT:=proc(f,s) local k;
t->evalf(Limit((-1)^k/k!*(k/t)^(k+1)*eval(diff(f,s$k),s=k/t),k=infinity)) 
end;
but later I corrected it so that it would work for a wider class of functions. Here is a simple example, 
inttrans[invlaplace](2/s^3,s,t);
                                  t^2
ILT(2/s^3,s)(10);
                             100.0000000
It can be plotted, 
plot(ILT(2/s^3,s),-1..1);
It can be also used with 3 arguments, where the 3rd argument is the number of digits. For example, 
ILT(2/s^10,s)(1000000);
                          .5511463845e49 
ILT(2/s^10,s,20)(1000000);
                     .55114638447971781305e49
Still, it works only for functions that can be differentiated quite simply. It doesn't work for most of test examples in Axel Vogt's post below. I'd like to thank Axel Vogt and Georgios Kokovidis for their comments below. Axel Vogt's program works for much wider class of functions than ILT. A good thing would be if a future version of Maple provided access to NAG Fortran library functions for inverse Laplace transform, C06LBF and C06LCF (Weeks' method) and an older procedure C06LAF (Crump's method). ______________ Alec Mihailovs http://mihailovs.com/Alec/
numeric inverse transform laplace

Please Wait...