Ask Your Question
2

How to find inverse laplace transform

asked 2017-02-18 13:13:43 -0600

ismon gravatar image

updated 2017-02-18 13:14:40 -0600

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 flag offensive close merge delete

2 answers

Sort by Β» oldest newest most voted
2

answered 2017-02-18 14:03:03 -0600

mforets gravatar image

updated 2017-02-22 03:44:17 -0600

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.

edit flag offensive delete link more

Comments

1

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

paulmasson gravatar imagepaulmasson ( 2017-02-19 17:51:44 -0600 )edit
1

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

kcrisman gravatar imagekcrisman ( 2017-02-20 08:17:22 -0600 )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?

mforets gravatar imagemforets ( 2017-02-22 03:18:13 -0600 )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.

kcrisman gravatar imagekcrisman ( 2017-02-22 18:28:09 -0600 )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.

mforets gravatar imagemforets ( 2017-02-23 05:04:57 -0600 )edit
0

answered 2017-02-19 06:28:58 -0600

ismon gravatar image

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?

edit flag offensive delete link more

Comments

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!

mforets gravatar imagemforets ( 2017-02-19 11:02:08 -0600 )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.

mforets gravatar imagemforets ( 2017-02-19 11:17:58 -0600 )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)').

mforets gravatar imagemforets ( 2017-02-22 03:06:50 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-02-18 13:13:43 -0600

Seen: 371 times

Last updated: Feb 22