Solving Laplace Inversion with Erfc

[matematikawan]

When trying to solve a pde using Laplace transform, I need to invert an expression of the form
[tex]\frac{\exp{(-as-b\sqrt{s})}}{s^2}[/tex]

A friend told me that Mathematica cannot invert such expression. I try using convolution but a bit loss when trying to evaluate the integral of Erfc(.).
[tex]L^{-1}\{\frac{\exp{(-as-b\sqrt{s})}}{s^2} \}[/tex]

[tex]=L^{-1} \{ \frac{e^{-as}}{s} \} * L^{-1} \{ \frac{e^{-b\sqrt{s}}}{s} \} [/tex]

[tex]=H(t-a) * Erfc(\frac{b}{2\sqrt{t}}\) [/tex]

[tex]=\int_0^t H(t-\tau-a) Erfc(\frac{b}{2\sqrt{\tau}}\) d\tau[/tex]

How do we proceed from here?

 

[matematikawan]

Is it all right if I proceed this way
[tex]=\int_{t-a}^t \mbox{Erfc}(\frac{b}{2\sqrt{\tau}}) d\tau[/tex]

[tex]= F(t) - F(t-a)[/tex] where F(x)=http://integrals.wolfram.com/index.jsp?expr=Erfc[b/2/Sqrt[x]]&random=false"


Any comment ?

 

[trambolin]

Only a suggestion but you can check the convergence of succesive Páde approximations for the exponential in your first post and maybe from analyticity you can have a conclusion.

 

[matematikawan]

Thanks trambolin for the suggestion on cross checking the answer. Although new to the method, I will give a try with this Páde approximations (some website write it as approximant).

It looks like there are scripts for mathematica at certain website that allow us to compute the coefficients of the Páde approximations. But a MATLAB script will greatly help because I'm more familier with this software. Does anyone seen such a MATLAB script before?

 

[trambolin]

actually if you type "doc pade" you can see what the function is, in MATLAB. But it works only for [itex]e^{sT}[/itex]. You can modify it a little bit further with some amount of boring work for your own function.

Actually the whole thing is nothing but the Taylor expansion and truncation. You can use even symbolic toolbox of MATLAB I guess.

 

[jackmell]

Is it all right if I proceed this way
[tex]=\int_{t-a}^t \mbox{Erfc}(\frac{b}{2\sqrt{\tau}}) d\tau[/tex]

[tex]= F(t) - F(t-a)[/tex] where F(x)=http://integrals.wolfram.com/index.jsp?expr=Erfc[b/2/Sqrt[x]]&random=false"Any comment ?

If t-a>0, then I believe it should be:

[itex]\int_{0}^{t-a} \text{erfc}\left(\frac{b}{2\sqrt{\tau}}\right)d\tau=F(t-a)-\lim_{\tau\to 0} F(\tau)[/itex]

since the antiderivative returned by Mathematica has a [itex]\sqrt{\tau}[/itex] in the denominator. Also, you can check all this numerically by just numerically integrating the inverse transform from say 1-k i to 1+k i for some reasonable range unless the integrand is badly-behaved. For example, when a=b=1 and t=5, just integrating it from 1-30 i to 1+30 i gives an answer in agreement to the antiderivative above correct to 5 decimal places.

 

[matematikawan]

Thanks everyone. My progress is a bit slow. However I do manage to write a script for computing the coefficients of the Pade approximant. The script gives the same result for exp(x) as found in the book Essentials of Pade Approximants by George A.Baker, Jr. Comments welcome.

function [a b] = pade(f,N);

% Compute the coefficients for an O(N)/O(N) Pade'
% approximant to the symbolic function f(x) about x=0
%
% Parameters:
% f - orginal function (using symbolic variable)
% N - highest order term in approximant
% y - point to expand about (number)
% a - coefficients of the numerator (decending)
% b - coefficients of the denominator (decending)

% Present version only expands around x=0
%
% Modified from: http://www.mathworks.com/matlabcentral/fileexchange/4388-pade-approximant"01/18/2004


Nc=2*N+1;

% Symbolic computation of first 2N+1 terms of the taylor series
fser=taylor(f,0,Nc);

% Extract the expansion coefficients
c=sym2poly(fser)';
fser=taylor(f,0,Nc);

% Extract the expansion coefficients
c=sym2poly(fser)';

% Reverse the order (ascending powers)
c=c(Nc:-1:1);

lhs=zeros(N);
for r=1:N
lhs(r,:)=c(r+1:N+r);
end

for r=1:N
rhs(r)=-c(N+r+1);
end

b=lhs\rhs';

b=b(N:-1:1);
b=[1; b]'; % b is denominator coefficients b1+b2*x+b3*x^2+...

a(1)=c(1);
for r=2:(N+1)
total=c(r);
for k=1:(r-1)
total=total+b(k+1)*c(r-k);
end
a(r)=total; % a is numerator coefficients
end