from sympy import *
import MATH280

x = Function("x")
t,s = symbols("t s",real=True,  nonnegative=true )

MATH280

Current sympy version 1.9 and hoping for a fix soon

A famously intense pen-and-paper technique from elementary ODEs class: Laplace Transforms. Sympy v1.9 has Laplace capabilities built-in, but with maddening shortcomings. This may be addressed in an upcomming new version. I see in the pre-release notes for sympy 1.10 "The Laplace transform is now rule-based, can transform a far wider range of expressions, and works better with undefined functions, including the ability to transform differential equations. "

Let's work with a 2nd order, linear, constant-coefficient, nonhomogeous IVP: $$ \ddot{x}+3\dot{x}+4x = \sin(2t), \;\;\;\;\;\; x(0)=1, \dot{x}(0)=0$$ 

ode = Eq(x(t).diff(t,2)+3*x(t).diff(t)+4*x(t) , sin(2*t) )
ode
$\displaystyle 4 x{\left(t \right)} + 3 \frac{d}{d t} x{\left(t \right)} + \frac{d^{2}}{d t^{2}} x{\left(t \right)} = \sin{\left(2 t \right)}$

laplace_transform() in sympy 1.9

First let's see how the built-in sympy laplace_transform() handles that equation:

laplace_transform(ode,t,s)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [3], in <module>
----> 1 laplace_transform(ode,t,s)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\sympy\integrals\transforms.py:1256, in laplace_transform(f, t, s, legacy_matrix, **hints)
   1253         else:
   1254             return type(f)(*f.shape, elements_trans)
-> 1256 return LaplaceTransform(f, t, s).doit(**hints)

File ~\AppData\Local\Programs\Python\Python310\lib\site-packages\sympy\integrals\transforms.py:180, in IntegralTransform.doit(self, **hints)
    174     raise IntegralTransformError(
    175         self.__class__._name, self.function, 'needeval')
    177 # TODO handle derivatives etc
    178 
    179 # pull out constant coefficients
--> 180 coeff, rest = fn.as_coeff_mul(self.function_variable)
    181 return coeff*self.__class__(*([Mul(*rest)] + list(self.args[1:])))

AttributeError: 'Equality' object has no attribute 'as_coeff_mul'

Look at that! 'Equality' object has no attribute. laplace_transform() would prefer that we assume the expression we enter is equal to 0:

Lx=laplace_transform(x(t).diff(t,2)+3*x(t).diff(t)+4*x(t) - sin(2*t),t,s)
Lx

(4*LaplaceTransform(x(t), t, s) + 3*LaplaceTransform(Derivative(x(t), t), t, s) + LaplaceTransform(Derivative(x(t), (t, 2)), t, s) - 2/(s**2 + 4),
 0,
 True)

We see the output above is a tuple that includes some conditional statements at the end. To hide that and see just the transform itself, use the option noconds=True

Lx=laplace_transform(x(t).diff(t,2)+3*x(t).diff(t)+4*x(t) - sin(2*t),t,s,noconds=True)
Lx
$\displaystyle 4 \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 3 \mathcal{L}_{t}\left[\frac{d}{d t} x{\left(t \right)}\right]\left(s\right) + \mathcal{L}_{t}\left[\frac{d^{2}}{d t^{2}} x{\left(t \right)}\right]\left(s\right) - \frac{2}{s^{2} + 4}$

Laplace Transform and Derivatives

Notice in that output another hassle: laplace_transform() ignores the single most useful property of the Laplace Transform: $$ {\cal L}\{x'(t)\} = s{\cal L}\{x(t)\} - x(0), \mbox{ and } {\cal L}\{x''(t)\} = s^2{\cal L}\{x(t)\} - sx(0)- x'(0).$$

laplace() in MATH280

We've reworked that and included an alternative laplace() in the MATH280 Module.

Solving an equation with Laplace Transforms in four steps:

1. take the transform of everything 

L = MATH280.laplace(ode,t,s)
L
$\displaystyle s^{2} \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 3 s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - s x{\left(0 \right)} + 4 \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - 3 x{\left(0 \right)} - \left. \frac{d}{d t} x{\left(t \right)} \right|_{\substack{ t=0 }} = \frac{2}{s^{2} + 4}$

2. plug in the initial conditions

The second step in solving the equation is to plug in the initial conditions:

L0=L.subs(x(0),1).subs(Subs(Derivative(x(t), t), t, 0),0)
L0
$\displaystyle s^{2} \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 3 s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - s + 4 \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - 3 = \frac{2}{s^{2} + 4}$

3. solve for the lapace transform of the solution function

The third step is to solve the resulting equation for the symbol ${\cal L}\{x(t)\}$

Lx=solve(L0,LaplaceTransform(x(t),t,s))
Lx

[(s**3 + 3*s**2 + 4*s + 14)/(s**4 + 3*s**3 + 8*s**2 + 12*s + 16)]

4. look up the laplace transform to determine the solution

The fourth and final step is to "look up" that complicated expression in the variable $s$ to determine the function of $t$ with that Laplace transform. The built-in sympy funtion inverse_laplace_transform() works fine for that, with the slight annoyance that it doesn't align perfectly with the list format of output from solve(). So I've included in MATH280 a little wrapper called laplaceInv() that extracts the zeroth entry from that list. I've noticed that this step can run really slowly.

sol = MATH280.laplaceInv(Lx,s,t)
sol
$\displaystyle - \frac{\cos{\left(2 t \right)} \theta\left(t\right)}{6} + \frac{\sqrt{7} e^{- \frac{3 t}{2}} \sin{\left(\frac{\sqrt{7} t}{2} \right)} \theta\left(t\right)}{2} + \frac{7 e^{- \frac{3 t}{2}} \cos{\left(\frac{\sqrt{7} t}{2} \right)} \theta\left(t\right)}{6}$

Something to notice is that every term as a $\theta(t)$. That's the unit step function: $\theta(t)=1$ if $t\gt 0$ and otherwise zero. Seems to be enforcing the assumption that our solution only makes sense for nonnegative time $t$.

We can plot the solution with sympy plot:

plot(sol,(t,0,20))

<sympy.plotting.plot.Plot at 0x2b67c954d90>

And just in case there was any doubt, we notice that dsolve() produces the same solution.

dsolve(ode,x(t),ics={x(0):1, x(t).diff(t).subs(t,0):0})
$\displaystyle x{\left(t \right)} = \left(\frac{\sqrt{7} \sin{\left(\frac{\sqrt{7} t}{2} \right)}}{2} + \frac{7 \cos{\left(\frac{\sqrt{7} t}{2} \right)}}{6}\right) e^{- \frac{3 t}{2}} - \frac{\cos{\left(2 t \right)}}{6}$

Discontinuous Forcing Functions

We haven't yet explored what I think is the reason to consider Laplace Transforms in the first place: discontinuous forcing functions.

A first-order model with a square wave switch function

Suppose we model a capacitor with an external voltage source that switches on at $t=2$ and off at $t=5$.
The equation, in some idealized units, could be $$\dot{x} - x = \theta(t-2) - \theta(t-5), \;\;\;\;\; x(0)=0 $$ where $\theta(t-a)$ is the Heaviside unit step function that turns on at $t=a$.

capmodel = Eq(x(t).diff(t)+x(t), Heaviside(t-2)-Heaviside(t-5))
capmodel
$\displaystyle x{\left(t \right)} + \frac{d}{d t} x{\left(t \right)} = - \theta\left(t - 5\right) + \theta\left(t - 2\right)$

We'll complete the four steps as in the previous example:

L=MATH280.laplace(capmodel,t,s)
L
$\displaystyle s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - x{\left(0 \right)} = \frac{e^{- 2 s}}{s} - \frac{e^{- 5 s}}{s}$ 
L0=L.subs(x(0),0)
L0
$\displaystyle s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) = \frac{e^{- 2 s}}{s} - \frac{e^{- 5 s}}{s}$ 
Ls=solve(L0,LaplaceTransform(x(t),t,s))
Ls

[(exp(3*s) - 1)*exp(-5*s)/(s*(s + 1))]
capsol=MATH280.laplaceInv(Ls,s,t)
capsol
$\displaystyle - \theta\left(t - 5\right) + \theta\left(t - 2\right) + e^{5} e^{- t} \theta\left(t - 5\right) - e^{2} e^{- t} \theta\left(t - 2\right)$

The functional form of the solution looks a little opaque, so lets make a plot of the solution together with the right hand side function:

plot(capsol,
     Heaviside(t-2)-Heaviside(t-5),
    (t,0,15),)

<sympy.plotting.plot.Plot at 0x2b67dc15330>

A second-order example with a delta function

Here's an example from the undergraduate ODE textbook by Blanchard, Devaney and Hall: $$\ddot{x} +2\dot{x}+3x= \delta(t-1)+3\delta(t-4), \;\;\;\;\; x(0)=0, \dot{x}(0)=0$$ 

ode2 = Eq(x(t).diff(t,2)+2*x(t).diff(t)+3*x(t),DiracDelta(t-1)+3*DiracDelta(t-4))
ode2
$\displaystyle 3 x{\left(t \right)} + 2 \frac{d}{d t} x{\left(t \right)} + \frac{d^{2}}{d t^{2}} x{\left(t \right)} = 3 \delta\left(t - 4\right) + \delta\left(t - 1\right)$ 
L=MATH280.laplace(ode2,t,s)
L
$\displaystyle s^{2} \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 2 s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - s x{\left(0 \right)} + 3 \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) - 2 x{\left(0 \right)} - \left. \frac{d}{d t} x{\left(t \right)} \right|_{\substack{ t=0 }} = e^{- s} + 3 e^{- 4 s}$ 
Ls=L.subs(x(0),0).subs(Subs(Derivative(x(t), t), t, 0),0)
Ls
$\displaystyle s^{2} \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 2 s \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) + 3 \mathcal{L}_{t}\left[x{\left(t \right)}\right]\left(s\right) = e^{- s} + 3 e^{- 4 s}$ 
Lx=solve(Ls,LaplaceTransform(x(t),t,s))
Lx

[(exp(3*s) + 3)*exp(-4*s)/(s**2 + 2*s + 3)]
sol2 = MATH280.laplaceInv(Lx,s,t)
sol2
$\displaystyle \frac{3 \sqrt{2} e^{4} e^{- t} \sin{\left(\sqrt{2} t \right)} \cos{\left(4 \sqrt{2} \right)} \theta\left(t - 4\right)}{2} + \frac{\sqrt{2} e e^{- t} \sin{\left(\sqrt{2} t - \sqrt{2} \right)} \theta\left(t - 1\right)}{2} - \frac{3 \sqrt{2} e^{4} e^{- t} \sin{\left(4 \sqrt{2} \right)} \cos{\left(\sqrt{2} t \right)} \theta\left(t - 4\right)}{2}$ 
plot(sol2,(t,0,20))

<sympy.plotting.plot.Plot at 0x2b67d9b30d0>

We can see in the plot above how the impulses at $t=1$ and $t=4$ energize the system.