Ask Your Question

Revision history [back]

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)$.

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)$. $5\cos(3(t-2))\mu(t-2)$.

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. Using the piecewise machinery, we can treat the result as an ordinary function. In the end of the day, it behaves as:

sage: InverseLaplace(5*s/(s^2 + 9), s, t) t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t) sage: plot(_, (t, -3, 3))


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

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

    Returns the list of summands of f.

    INPUT: 

    * "f" - symbolic expression

    EXAMPLES:

    Mononmials 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can treat the result as an ordinary function. In the end of the day, it behaves as:

 sage: InverseLaplace(5*s/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))


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

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

    Returns the list of summands of f.

    INPUT: 

    * "f" - symbolic expression

    EXAMPLES:

    Mononmials 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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 inverse_laplace(5*s*exp(-2*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. Using the piecewise machinery, we can treat the result as an ordinary function. In the end of the day, it behaves as:

sage: InverseLaplace(5*s/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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

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

    Returns the list of summands of f.

    INPUT: 

    * "f" - symbolic expression

    EXAMPLES:

    Mononmials 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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*exp(-2*s)/(s^2 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. Using the piecewise machinery, we can treat the result as an ordinary function. In the end of the day, it behaves as:

sage: InverseLaplace(5*s/(s^2 InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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

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

    Returns the list of summands of f.

    INPUT: 

    * "f" - symbolic expression

    EXAMPLES:

    Mononmials 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can treat the result as an ordinary function. In the end of the day, it behaves as:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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

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

    Returns the list of summands of f.

    INPUT: 

    * "f" - symbolic expression

    EXAMPLES:

    Mononmials 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can treat the result as produce an ordinary function. function as a result. In the end of the day, it behaves as:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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

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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves as:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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

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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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,
    where: 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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves as:like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

First we need this one, to treat cases involving a sum in the numerator (not in the OP's question):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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

Finally, there is a comparison with the 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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 in time and also in the output result (this was already reported here):

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

while:

%%time

g = 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

g

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))

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

Finally, there is a comparison with the 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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 in time and also in the output result (this was already reported here):

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

while:

%%time

g 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

g
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))

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

Finally, there is a comparison with the 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 piecewise construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))

        return sum(op_ilt)

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):

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

while:

%%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))

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. Using the piecewise machinery, we can produce an ordinary function as a result. In the end of the day, it behaves like:

sage: InverseLaplace(5*s*exp(-2*s)/(s^2 + 9), s, t)
t |--> 5*cos(3*t - 6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)
sage: plot(_, (t, -3, 3))

$$ \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 piecewiseheaviside construct.

def heaviside_step_function(t, a=0):
    """Heaviside step function shifted to a. 

    Implements the piecewise function: mu(t) = 1 if t < a and 0 if t >= a.

    INPUT: 

    * "t" - variable of the output function.

    * "a" - (default = 0), time-shift mu(t) -> mu(t-a) of the canonical Heaviside function mu(t). 

    EXAMPLES: 

    sage: var('t')
    sage: plot(heaviside_step_function(t), (t, -3, 3)) + plot(heaviside_step_function(t, 2), (t, -3, 3), color='red')

    Negative values of a are also allowed:

    sage: plot(heaviside_step_function(t, -1), (t, -2, 2)) 

    NOTES: 

    - The extra interval [a, a] is there to define the function at zero. Don't know if there is a way to 
    define half-closed intervals.

    - An improvement (for readability when computing with Laplace transforms) would be to wrap this into a 
    function \mu(t), so that it works like:

    sage: var('t')
    sage: heaviside_step_function(t, 2);
    mu(t-2)
    sage: _(t=1), _(t=2)
    (0, 1)
    """

    # instantiate PW function
    mu(t) = piecewise([[(-oo, a), 0], [[a, a], 1], [(a, oo), 1]], var=t)

    return mu
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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); 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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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_step_function(t, -a))
op_ilt.append(fi(t=t+a)*heaviside(t+a))

        return sum(op_ilt)

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):

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

while:

%%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))

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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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 sum(op_ilt)
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):

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

while:

%%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))

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)
0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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

while:

%%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))

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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d = op.match(w1*exp(w0*s))

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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

while:

%%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))

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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d0 = op.match(exp(w0*s))
            d = op.match(w1*exp(w0*s))
op.match(w1*exp(w0*s)) if d0 == None else d0

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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

while:

%%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))

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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d0 = op.match(exp(w0*s))
            d = op.match(w1*exp(w0*s)) if d0 == None else d0

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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

while: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). $$

While

%%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))

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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d0 = op.match(exp(w0*s))
            d = op.match(w1*exp(w0*s)) if d0 == None else d0

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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). $$

While 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))

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) 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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d0 = op.match(exp(w0*s))
            d = op.match(w1*exp(w0*s)) if d0 == None else d0

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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):

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))

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, 0, 10)

    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

    """

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

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

    # check if time-shifting applies
    got_time_shifted = True if num.has(exp(w0*s)) and not den.has(exp(w0*s)) 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:

            d0 = op.match(exp(w0*s))
            d = op.match(w1*exp(w0*s)) if d0 == None else d0

            if d == None:

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

            elif w0 in d:

                a = d[w0]

                # 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))

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, 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)
6)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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)
1))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    sage: plot(h, 0, 10)
-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).expand() 
    den = F.denominator()(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 num.has(exp(w0*s)) got_time_shifted_num and not den.has(exp(w0*s)) 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:

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

            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))

            elif else:

                a = d[w0] if w0 in d:

                a = d[w0]
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)        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))

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); t)

    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 |--> (2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); t)

    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).expand() F.numerator()(s=s).expand() 
    den = F.denominator()(s).expand()
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))

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)*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 2), t|-->1 on {2}, 
    t|-->1 on (2, +oo); 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))*piecewise(t|-->0 on (-oo, 1), t|-->1 on {1}, t|-->1 on (1, +oo); 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))

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.

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 I doesn't seem obvious how to bring it back to Sage via SR(f) for evaluation, plotting, etc.

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.

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 I it doesn't seem obvious how to bring it back to Sage via for evaluation, plotting, etc, as SR(f) for evaluation, plotting, etc.

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.

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 (For the record, this example is run in a Jupyter notebook with SageMath 7.5.1 kernel. 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.