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()
if 'add' not in str(f.operator()):
# no addition/substraction found
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
.