Revision history [back]
 | 1 | initial version |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(t−2).
| | 2 | No.2 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(t−2). 5cos(3(t−2))μ(t−2).
| | 3 | No.3 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 4 | No.4 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 5 | No.5 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 6 | No.6 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 7 | No.7 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 8 | No.8 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 9 | No.9 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 10 | No.10 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 11 | No.11 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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)| | 12 | No.12 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 13 | No.13 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 14 | No.14 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(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 muWe 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 15 | No.15 Revision |
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, L−1(e−asF(s))=f(t−a)μ(t−a) (where μ(t) is the Heaviside step function), to conclude that the answer is 5cos(3(t−2))μ(t−2).
For exercising, below is an attempt to mechanize this rule. Using the In the end of the day, it behaves like:piecewise machinery, we can produce an ordinary function as a result.
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 piecewiseheavisideconstruct.
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 muWe 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 16 | No.16 Revision |
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)
fKicking 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 17 | No.17 Revision |
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 fKicking 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 18 | No.18 Revision |
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 fKicking 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 19 | No.19 Revision |
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 fKicking 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 loopCompare 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 sThe 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 loopwhile:
%%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 sand
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))| | 20 | No.20 Revision |
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 fKicking 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 loopCompare 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 sThe 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 loopwhile: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 sand
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))| | 21 | No.21 Revision |
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 fKicking 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 loopCompare 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 sThe 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 loopand 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 sand
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))| | 22 | No.22 Revision |
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))
| | 23 | No.23 Revision |
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))
| | 24 | No.24 Revision |
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))
| | 25 | No.25 Revision |
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))
| | 26 | No.26 Revision |
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))
| | 27 | No.27 Revision |
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.
| | 28 | No.28 Revision |
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.
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. | | 29 | No.29 Revision |
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.
