How to find inverse laplace transform

edit

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.