# Problem with Inverse Laplace Transform in Mathematica

I am trying to calculate the following inverse Laplace transform

$$f(t) = \mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\left( \frac{1}{s+b}\right)^{n} \right](t)$$

with $$s$$ being the variable, $$a,b >0$$ and $$n\in N$$. I know I can use the properties of the inverse transform to calculate this from the convolution

$$\mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\right]*\mathcal{L}^{-1}\left[\left( \frac{1}{s+b}\right)^{n} \right](t)$$

and the fact that:

$$\mathcal{L}^{-1}\left[ \left(\frac{1}{s+a}\right)^{n+1}\right] = \frac{t^n e^{-at}}{n!}\Theta(t)$$

Combining this together this gives for the final result of the inverse transform:

$$f(t) = \frac{\sqrt{\pi t (b-a)}}{2 \Gamma(n+1)} e^{-\frac{1}{2}(a+b)t} \left[I_{n-1/2}\left( \frac{1}{2}(b - a) t\right)+I_{n+1/2}\left( \frac{1}{2} (b - a) t \right) \right] ,$$

where $$I_{\nu}(x)$$ is the Modified Bessel function of the first kind. However, if I try to check whether or not this is correct by taking the Inverse transform I get an expression that only gives the initial form if $$b>a$$, but not otherwise

I cannot understand what the problem is, or if I am missing something important.

Here is the Mathematica code I am using:

fun1 = 1/(s + a)^(NN + 1)

fun2 = 1/(s + b)^NN

fun = fun1*fun2

baseTransform = 1/(s + x)^M

baseInverse = InverseLaplaceTransform[baseTransform, s, t1]

fN = baseInverse /. {M -> (NN + 1), x -> a, t1 -> t - tau}

gN = baseInverse /. {M -> NN, x -> b, t1 -> tau}

sol = Refine[Integrate[fN*gN, {tau, 0, t}], {NN \[Epsilon] Integer, NN > 0, t \[Epsilon] Real, t > 0}]

funCheck1 = Simplify[LaplaceTransform[sol, t, s],Assumptions -> {a < b}]

funCheck2 = Simplify[LaplaceTransform[sol, t, s],Assumptions -> {a > b}]

• Welcome to the Mathematica Stack Exchange :) Please add the Wolfram Language code you are using to the question - typing out those equations is a lot of work! – Carl Lange Mar 24 '20 at 11:29
• @Carl LangeL: Done, sorry for not doing before, the code is now added at the end of the question – user2918065 Mar 24 '20 at 12:54

If we use:

L1 = InverseLaplaceTransform[1/(s + a)^(n + 1), s, t];
L2 = InverseLaplaceTransform[1/(s + b)^n, s, t];
L12 = FullSimplify[Integrate[(L1 /. t -> (t - \[Tau]))*(L2 /. t -> \[Tau]), {\[Tau], 0,t}, Assumptions -> n > 0], t > 0]

f[a_, b_, n_, t_] := InverseLaplaceTransform[1/(s + a)^(n + 1) 1/(s + b)^n, s, t] // N;
g[a_, b_, n_, t_] := ((-a + b)^(1/2 - n) E^(-(1/2) (a + b) t) Sqrt[\[Pi]] t^(
1/2 + n) (BesselI[-(1/2) + n, 1/2 (-a + b) t] + BesselI[1/2 + n, 1/2 (-a + b) t]))/(2 Gamma[1 + n]) // N


Then we can check and works fine for a > b and b < a:

(Table[g[a, b, n, 1] - f[a, b, n, 1] /. Indeterminate -> 0, {a, -5, 5}, {b, -5, 5}, {n, -5, 5}] // Quiet // Chop) // MatrixForm(* for t = 1 *)

(*ALL Zeros *)

• Thanks for the answer but somehow this is not what I am looking for. Numerically the solution looks fine, but if you try to calculate it symbolically you will see that the inverse of f is not L12 if a>b. I cannot understand why is that so and if there is some hidden problem I cannot understand. One more quick thing, you should change your answer "Then we can check and works fine for a > b and b < a" should read: "Then we can check and works fine for a > b and b > a" (or it would be the same case). – user2918065 Mar 24 '20 at 17:57
• For me formula works fine for: a > b and b > a. If you don't like my answer, wait for a better answer. :P – Mariusz Iwaniuk Mar 24 '20 at 18:04
• when you say your formula "works" it is a bit ambiguous. Numerically, I agree with you and it is easy to test. However, why if I then use LaplaceTransform on the solution Mathematica does not provide the initial value of the transform for a>b? Is it possibly a bug, or is there something important missing? – user2918065 Mar 25 '20 at 19:56