Search code examples
pythontransformsympyintegral

Laplace transform of unknown integral (function of time)


I need to calculate the Laplace transform of an integral function. It seems that sympy is not yet able to understand that.

Assuming the following:

from sympy import *
s, t = symbols('s t')
I = Function('I')(t)
eq1 = integrate(I, t)
transforms.laplace_transform(eq1, t, s)

The solution should be: I(s) / s

However, sympy gives: LaplaceTransform(Integral(I(t), t), t, s)

It seems to be an open issue Issue 7219. Is there any work around?


Solution

  • It seems that the issue hasn't been fixed yet.

    However, we can give a "crappy workaround" based on the "crappy implementation" of Eric Wieser given for derivatives. Note, however, that the original snippet doesn't seem to work for derivatives either, because the internal representation of higher order derivatives seems to have changed since the snippet was posted.

    Here's my "crappy" workaround that catches only the simplest of cases (derivatives only with respect to t, indefinite integrals only with respect to t, where t is the variable on which the Laplace transform acts):

    from sympy import *
    
    def laplace(e, t, s):
        """Hacked-up Laplace transform that handles derivatives and integrals
    
        Updated generalization of https://github.com/sympy/sympy/issues/7219#issuecomment-154768904
        """
    
        res = laplace_transform(e, t, s, noconds=True)
        wf = Wild('f')
        lw = LaplaceTransform(wf, t, s)
    
        for exp in res.find(lw):
            e = exp.match(lw)[wf]
            args = e.args
    
            if isinstance(e, Derivative):
                # for derivative check that there's only d/dt^n with n>0
                if len(args) == 2 and args[1][0] == t:
                    n = args[1][1]
                    if n > 0:
                        newexp = s**n * LaplaceTransform(e.args[0], t, s)
                    res = res.replace(exp, newexp)
    
            elif isinstance(e, Integral):
                # for integral check that there's only n consecutive indefinite integrals w.r.t. t
                if all(len(arg) == 1 and arg[0] == t for arg in args[1:]):
                    newexp = s**(-len(args[1:])) * LaplaceTransform(args[0], t, s)
                    res = res.replace(exp, newexp)
    
            # otherwise don't do anything
    
        return res
    
    x = Function('x')
    s,t = symbols('s t')
    print(laplace(Derivative(x(t), t, 3), t, s))
    print(laplace(Integral(Integral(x(t), t), t), t, s))
    

    The above outputs

    s**3*LaplaceTransform(x(t), t, s)
    LaplaceTransform(x(t), t, s)/s**2
    

    as expected. Using your specific example:

    I = Function('I')(t)
    eq1 = integrate(I, t)
    LI = laplace(eq1, t, s)
    print(LI)
    

    we get

    LaplaceTransform(I(t), t, s)/s
    

    which is the correct representation of "I(s)/s" that you expected.


    The way the above workaround works is that it matches the arguments of the LaplaceTransform and checks if there's a pure Derivative or Integral inside. For Derivative we check that there's only differentiation with respect to t; this is what Eric's original workaround did, but while his code seems to have expected args of the form Derivative(x(t), t, t, t), the current representation of derivatives is Derivative(x(t), (t,3)). This is why handling this use case had to be changed.

    As for Integrals, the representation is similar to the original one: Integral(x(t), t, t) is a double integral. I still had to tweak Eric's original, because the args of this expression contain tuples for each integral rather than a scalar t, in order to accommodate definite integrals. Since we only want to handle the no-brainer case of indefinite integrals, I made sure that there's only indefinite integrals and only with respect to t.

    If the argument of the LaplaceTransform is anything else, the expression is left alone.