Using Laplace Transforms to Solve a Linear Differential Equation in SymPy

Consider solving

\begin{equation*} 3\ddot{x}+30\dot{x}+63x=4\dot{g}(t)+6g(t) \end{equation*}

in Jupyter where

\begin{equation*} g(t)=u_s(t) \end{equation*}

is the unit step function (Heaviside Function) and \(x(0) = 4\) and \(\dot{x}(0)=7\).

%matplotlib inline
 from sympy.integrals.transforms import laplace_transform
 from sympy.integrals.transforms import inverse_laplace_transform
 from sympy import *
 import sympy as sympy
 from sympy.abc import a, t, x, s, X, g, G
 init_printing(use_unicode=True)

x, g, X = symbols('x g X', cls = Function)

x0 = 4
 v0 = 7

g = Heaviside(t) #This is just the step function

laplace_transform(Heaviside(t), t, s) # Checking the performance of the Laplace Transform Module

(1/s, 0, True)

laplace_transform(DiracDelta(t), t, s) # Checking the performance... this is wrong. Should get 1.

(1/2, -oo, True)

laplace_transform(DiracDelta(t - 2), t, s) # right result

(exp(-2*s), -oo, True)

laplace_transform(DiracDelta(t - a), t, s)[0] # right result (theta(t) is a Heaviside function representation in SymPy)

\begin{equation*} \left(- \theta\left(a e^{i \pi}\right) + 1\right) e^{- a s} \end{equation*}

laplace_transform(DiracDelta(t - a), t, s)[0].subs(a,2) # Looks good

\begin{equation*} e^{- 2 s} \end{equation*}

laplace_transform(DiracDelta(t - a), t, s)[0].subs(a,1) # Looks good

\begin{equation*} e^{- s} \end{equation*}

laplace_transform(DiracDelta(t - a), t, s)[0].subs(a,0.0000001) # Looks to converge

\begin{equation*} e^{- 1.0 \cdot 10^{-7} s} \end{equation*}

laplace_transform(DiracDelta(t - a), t, s)[0].subs(a,-0.0000001) #So... here is the illustration.

\begin{equation*} 0 \end{equation*}

What's happening here is that SymPy currently takes the position that half the Dirac delta happens before zero, half after, so the result should only be half as big. I think I'm in a rather large crowd in saying that this isn't proper. It makes some theoretical sense, and is a wonderful math debate. However, given convention says that \(\delta(t)\) is fully captured by a Laplace transform with a result of \(1\) (Mathematica, Maple, Matlab, every System Dynamics, Controls, and Signal Processing book I've ever read), SymPy is practically wrong. I'm hoping that they will change their minds. I am now a bit skeptical about using SymPy for my math work as the results of a simple conventional application don't match expected behavior consistent with other codes.[1]_

So we note that SymPy isn't taking the Laplace Transform properly here, so we need to avoid using this result. (we should have gotten 1) Valid as of 0.7.6.1

print(sympy.__version__)

0.7.6.1

The first line below would work if SymPy performed the Laplace Transform of the Dirac Delta correctly. Short of that, we manually insert the Laplace Transform of \(g(t)\) and \(\dot{g}(t)\) where \(g(t)=u(t)\).

#eom  = Eq(3*(s**2 * X(s)-s*x0-v0)+30 * (s*X(s)-x0)+ 63 * X(s), laplace_transform( 4 * diff(g,t) + 6 * g, t, s, noconds = True))

 eom  = Eq(3*(s**2 * X(s)-s*x0-v0)+30 * (s*X(s)-x0)+ 63 * X(s), 4 * 1 + 6 * 1/s)
 eom

\begin{equation*} 3 s^{2} X{\left (s \right )} + 30 s X{\left (s \right )} - 12 s + 63 X{\left (s \right )} - 141 = 4 + \frac{6}{s} \end{equation*}

Xofs = solve(eom,X(s))
Xofs[0]

\begin{equation*} \frac{12 s^{2} + 145 s + 6}{3 s \left(s^{2} + 10 s + 21\right)} \end{equation*}

soln = inverse_laplace_transform(Xofs[0],s,t)
soln

\begin{equation*} \frac{\theta\left(t\right)}{84} \left(8 e^{7 t} + 749 e^{4 t} - 421\right) e^{- 7 t} \end{equation*}

soln_simp = expand(soln)
soln_simp

\begin{equation*} \frac{2 \theta\left(t\right)}{21} + \frac{107 \theta\left(t\right)}{12} e^{- 3 t} - \frac{421 \theta\left(t\right)}{84} e^{- 7 t} \end{equation*}

N(soln_simp,5)

\begin{equation*} 0.095238 \theta\left(t\right) + 8.9167 e^{- 3 t} \theta\left(t\right) - 5.0119 e^{- 7 t} \theta\left(t\right) \end{equation*}

Note that \(\theta(t)\) is SymPy's notation for a step function. This simply means the answer can't be used before \(t=0\).

[1]	This has been fixed in later revisions.

Comments