# How to find inverse laplace transform

Hi!

I'm trying to find the inverse laplace transform from the following equation: (5 * s * e^(-2*s))/(s^2 + 9).inverse_laplace(s,t),

although sage responds with:

ilt((5 * s * e^(-2*s))/(s^2 + 9)).

How do I find the inverse laplace transform from the equation above?

Sincerly Simon

edit retag close merge delete

Sort by » oldest newest most voted

In the docstring (inverse_laplace?), we learn that ilt is returned when no explicit inverse Laplace transform is found. Let's try to simplify it a bit:

sage: inverse_laplace(5*s/(s^2 + 9), s, t)
5*cos(3*t)


We can apply the time shifting property, $\mathcal{L}^{-1}(e^{-as}F(s)) = f(t-a)\mu(t-a)$ (where $\mu(t)$ is the Heaviside step function), to conclude that the answer is $5\cos(3(t-2))\mu(t-2)$.

For exercising, below is an attempt to mechanize this rule. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)


$$\newcommand{\Bold}[1]{\mathbf{#1}}5 \, \cos\left(3 \, t - 6\right) H\left(t - 2\right).$$

Finally, there is a comparison with SymPy's homonymous inverse_laplace_transform.

First we need this one, to treat cases involving a sum in the numerator (not in the OP's question though):

def break_into_summands(f):
"""Break into summands a symbolic expression.

Returns the list of summands of f.

INPUT:

* "f" - symbolic expression

EXAMPLES:

Monomials are transformed into lists with one element:

sage: break_into_summands(x*2^x)
[2^x*x]

If the expression contains terms that are factored, tries to multiply them out:

sage: f = (5*x*exp(-2*x) - 25*x^4)*sin(2*pi*x) - 1/x
sage: break_into_summands(f)
[-25*x^4*sin(2*pi*x), 5*x*e^(-2*x)*sin(2*pi*x), -1/x]

"""

f = f.expand()

return [f]
else:
return f.operands()


To define the Heaviside step function, we use the built-in heaviside construct.

We are ready to define the extension:

def InverseLaplace(F, s, t):
"""Inverse Laplace transform, with time-shifting property.

Attempts to obtain an explicit inverse Laplace transform of F(s) by using the time-shifting property,
F(s) := G(s)*exp(-as) -> g(t-a)*mu(t-a),
with mu(t) being the Heaviside step function, and a is real.

It is assumed that F(s) = num(s)/den(s) is a strictly proper real-rational function with any number of time-shifts,
namely: num(s) = \sum_{k} n_k s^k exp(-a_k*s) and den(s) = \sum_{k} d_k s^k.

See inverse_laplace for further details.

EXAMPLES:

sage: var('t s')
sage: F(s) = (5*s*exp(-2*s))/(s^2+9)
sage: f = InverseLaplace(F, s, t); f

t |--> 5*cos(3*t - 6)*heaviside(t - 2)

sage: plot(f, -3, 3)

More than one real exponential in the numerator can be handled:

sage: H(s) = (3*s^2*exp(-2*s) - exp(-s))/(s^3+1)
sage: h = InverseLaplace(H, s, t); h

t |--> 1/3*((sqrt(3)*sin(1/2*sqrt(3)*(t - 1)) - cos(1/2*sqrt(3)*(t - 1)))*e^(1/2*t - 1/2) + e^(-t + 1))*heaviside(t - 1) + (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*heaviside(t - 2)

sage: plot(h, -3, 3)

NOTES:

- We slightly change the convention using exp(as), because has(exp(-w0*s)) doesn't work. For instance,
sage: expr = exp(-2*s)*5*s; w0 = SR.wild(0)
sage: expr.has(exp(w0*s))
True

but

sage: expr.has(exp(-w0*s))
False

- There is an awkward chain of "matches", again related to the behavior of wild cards and the identity. Don't know if there is a better way.

"""

w0 = SR.wild(0); w1 = SR.wild(1);

num = F.numerator()(s=s).expand()
den = F.denominator()(s=s).expand()

# check if time-shifting applies
got_time_shifted_num = (num.has(exp(w0*s)) or num.has(exp(s)))
got_time_shifted_den = (den.has(exp(w0*s)) or den.has(exp(s)))

got_time_shifted = True if got_time_shifted_num and not got_time_shifted_den else False

if not got_time_shifted:

return inverse_laplace(F, s, t)

else:

# decompose into summands (+, -)
num_operands = break_into_summands(num)

# collect inverse Laplace transform of each summand
op_ilt = []

for op in num_operands:

# this is because wild cards don't seem to like multiplication by 1 (see note above)
d = op.match(exp(w0*s))
if d == None:
d = op.match(exp(s))
if d == None:
d = op.match(w1*exp(s))
if d == None:
d = op.match(w1*exp(w0*s))

if d == None:

op_ilt.append(inverse_laplace(op/den, s, t))

else:

a = d[w0] if w0 in d else 1

# compute ILT of the factor in front of the exponential
fi = inverse_laplace(op.content(exp(a*s))/den, s, t)

# apply time shift and multiply by Heaviside step function
op_ilt.append(fi(t=t+a)*heaviside(t+a))

f(t) = sum(op_ilt)

return f


### Kicking the tires

Let us see if we can do this more easily with SymPy.

First, the OP's function using InverseLaplace:

timeit('InverseLaplace((5*s*exp(-2*s))/(s^2+9), s, t)')
125 loops, best of 3: 4.65 ms per loop


Compare to:

%%time

from sympy.integrals.transforms import inverse_laplace_transform
from sympy import exp, Symbol
from sympy.abc import s, t

f = inverse_laplace_transform((5*s*exp(-2*s))/(s^2+9), s, t);

CPU times: user 1.14 s, sys: 5.92 ms, total: 1.15 s
Wall time: 1.14 s


The result takes a bit more computation time but is nicely wrapped into the Heaviside function:

f
5*cos(3*t - 6)*Heaviside(t - 2)

# plot is also easy:
plot(SR(f), (t, 0, 3))


Let's say this is a draw.

In this second example, inverse_laplace_transform behaves worse also in the output result (this was already reported here). First, observe that:

var('s t')
H(s) = (3*s^2*exp(-2*s) - exp(-s))/(s^3+1)
timeit('InverseLaplace(H, s, t)')
25 loops, best of 3: 12.4 ms per loop


and produces

$$\newcommand{\Bold}[1]{\mathbf{#1}}t \ {\mapsto}\ \frac{1}{3} \, {\left({\left(\sqrt{3} \sin\left(\frac{1}{2} \, \sqrt{3} {\left(t - 1\right)}\right) - \cos\left(\frac{1}{2} \, \sqrt{3} {\left(t - 1\right)}\right)\right)} e^{\left(\frac{1}{2} \, t - \frac{1}{2}\right)} + e^{\left(-t + 1\right)}\right)} H\left(t - 1\right) + {\left(2 \, \cos\left(\frac{1}{2} \, \sqrt{3} {\left(t - 2\right)}\right) e^{\left(\frac{1}{2} \, t - 1\right)} + e^{\left(-t + 2\right)}\right)} H\left(t - 2\right).$$

On the other hand,

%%time

h = inverse_laplace_transform((3*s^2*exp(-2*s) - exp(-s))/(s^3+1), s, t);

CPU times: user 35.7 s, sys: 106 ms, total: 35.8 s
Wall time: 35.8 s


and

h

4*(I*(-3*sqrt(3)*exp(3*t/2 + 1/2)*sin(sqrt(3)*t/2 - sqrt(3) + pi/6)*Heaviside(t - 2)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2) + 3*sqrt(3)*exp(3*t/2 + 1/2)*sin(sqrt(3)*t/2 - sqrt(3) + pi/6)*Heaviside(t - 2)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2) - sqrt(3)*exp(3*t/2 + 1)*cos(sqrt(3)*(t - 1)/2)*Heaviside(t - 1)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2) + sqrt(3)*exp(3*t/2 + 1)*cos(sqrt(3)*(t - 1)/2)*Heaviside(t - 1)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2)) + 3*sqrt(3)*exp(3*t/2 + 1/2)*sin(-sqrt(3)*t/2 + pi/3 + sqrt(3))*Heaviside(t - 2)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2) + 3*sqrt(3)*exp(3*t/2 + 1/2)*sin(-sqrt(3)*t/2 + pi/3 + sqrt(3))*Heaviside(t - 2)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2) - sqrt(3)*exp(3*t/2 + 1)*sin(sqrt(3)*(t - 1)/2)*Heaviside(t - 1)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2) - sqrt(3)*exp(3*t/2 + 1)*sin(sqrt(3)*(t - 1)/2)*Heaviside(t - 1)*gamma(-1/2 - sqrt(3)*I/2)*gamma(-1/2 + sqrt(3)*I/2)*gamma(3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2) + 9*exp(7/2)*Heaviside(t - 2)*gamma(-3/2 - sqrt(3)*I/2)*gamma(-3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2) - 3*exp(5/2)*Heaviside(t - 1)*gamma(-3/2 - sqrt(3)*I/2)*gamma(-3/2 + sqrt(3)*I/2)*gamma(5/2 - sqrt(3)*I/2)*gamma(5/2 + sqrt(3)*I/2))*exp(-t - 3/2)*cosh(sqrt(3)*pi/2)**2/(3*pi**2*(3 - sqrt(3)*I)*(3 + sqrt(3)*I))


#### Using Giac/XCAS

This requires that you have installed Giac/XCAS (and that it is accessible through the command line). (For the record, this example is run in a Jupyter notebook with SageMath 7.5.1 kernel, Mac OSX 10.11.6, 2,2 GHz Intel Core i7).

Let's go straight to the 'hard' example:

%%time

# with giac
f = giac('invlaplace((3*s^2*exp(-2*s) - exp(-s))/(s^3+1), s, t)')
f

CPU times: user 714 µs, sys: 812 µs, total: 1.53 ms
Wall time: 3.09 ms


and we get:

$$\newcommand{\Bold}[1]{\mathbf{#1}}"(e^{-t+2}+2 \cos\left(\frac{\sqrt{3}}{2} \cdot (t-2)\right) e^{\frac{1}{2} \cdot (t-2)}) \mathrm{Heaviside}\left(t-2\right)+(-\frac{1}{3} \cdot e^{-t+1}+\frac{1}{3} \cdot \cos\left(\frac{\sqrt{3}}{2} \cdot (t-1)\right) e^{\frac{1}{2} \cdot (t-1)}-\frac{1}{\sqrt{3}} \cdot \sin\left(\frac{\sqrt{3}}{2} \cdot (t-1)\right) e^{\frac{1}{2} \cdot (t-1)}) \mathrm{Heaviside}\left(t-1\right)"$$

Executed in a loop, we get 125 loops, best of 3: 2.98 ms per loop, hence 4x faster than the Sage implementation, and more than 12000x faster than the SymPy call.

Here type(f) = <𝚌𝚕𝚊𝚜𝚜'𝚜𝚊𝚐𝚎.𝚒𝚗𝚝𝚎𝚛𝚏𝚊𝚌𝚎𝚜.𝚐𝚒𝚊𝚌.𝙶𝚒𝚊𝚌𝙴𝚕𝚎𝚖𝚎𝚗𝚝'>, but unfortunately it doesn't seem obvious how to bring it back to Sage for evaluation, plotting, etc, as SR(f) returns TypeError: unable to convert (exp(-t+2)+2*cos(sqrt(3)/2*(t-2))*exp(1/2*(t-2)))*Heaviside(t-2)+(-1/3*exp(-t+1)+1/3*cos(sqrt(3)/2*(t-1))*exp(1/2*(t-1))-1/sqrt(3)*sin(sqrt(3)/2*(t-1))*exp(1/2*(t-1)))*Heaviside(t-1) to a symbolic expression.

more

1

@mforets if you incorporate this into the Sage source code, the cc me and I'll review it!

( 2017-02-19 17:51:44 -0500 )edit
1

Yes, please do try to get that in if it's reasonably robust (evaluation etc.).

( 2017-02-20 08:17:22 -0500 )edit

@paulmasson and @kcrisman I appreciate your feedback! I'll try to make it my 1st ticket ;) .. but as it is above it's not great, for instance with this one: $exp(-s+1)/s$. Also for consistency we need to address the direct Laplace with heaviside/time-shift; this feature was already requested here. A higher-level question/comment: as you know this is already solved in other systems that Sage can interface with, not only in SymPy as shown above, I'm thinking in Giac/XCAS, which is fast and gives 'better' answers in this particular setting as some experiments show. So what about calling Giac/XCAS whenever the algorithm in Maxima returns no explicit inverse Laplace transform?

( 2017-02-22 03:18:13 -0500 )edit

Interesting. But unfortunately most people won't have that installed (I don't, for instance). At least adding this to documentation would be good; I'm not sure what to do long-term with non-standard packages for such behavior.

( 2017-02-22 18:28:09 -0500 )edit
1

ok thanks, that's what I thought, but actually I'm confused since there is a relatively recent closed ticket make giac a standard package ? I've created ticket #22422.

( 2017-02-23 05:04:57 -0500 )edit

Yes, I see your point. Thank you.

Although, I hoped that I shouldn't rely on analysis, since I would like to use sage to calculate the answer to compare with my own calculations and analysis. With MATLAB for instance, the answer can be given without doing any analysis.

Is this possible with sage?

more

well, a fairly direct approach is to feed inverse_laplace with simpler expressions and from the "exterior" apply systematically the target property that we are interested in. see my edited answer above, comments welcome!

( 2017-02-19 11:02:08 -0500 )edit

... as it is often the case, it is also useful to check if SymPy provides the functionality that you are looking for; in this case yes (see the doc). I edited my answer with a comment on this.

( 2017-02-19 11:17:58 -0500 )edit
1

yet another way of doing this in sage, it requires that you have installed the system Giac/XCAS: sage: giac('invlaplace((5*s*exp(-2*s))/(s^2+9), s, t)').

( 2017-02-22 03:06:50 -0500 )edit