Using Laplace Transforms in Python to Solve Differential Equations

Photo by Jeremy Bishop on Unsplash

Last time, we played with the Laplace transform and its inverse in Python. Today, we’ll use that knowledge to solve differential equations.

Suppose we have the differential equation

for the unknown function 𝑦(𝑡). This equation describes a forced oscillator with friction in physics. As initial conditions, we’ll choose 𝑦(0)=𝑦′(0)=0.

The Laplace transform offers a most convenient method to solve this kind of equation. First of all, look what happens, if we Laplace transform the second derivative of our unknown function:

where we have used partial integration from the first to the second line. So we can replace the second derivative by multiplying with p and subtracting with the initial condition for the first derivative. For 𝐿(𝑦′) we just do the same and get:

So we get

This allows us to Laplace transform the whole differential equation. Let’s switch to Python and start a Jupyter notebook. Define the symbols and the differential equation, along with the unevaluated Laplace transform:

from sympy import *

t, p = symbols('t, p')
y = Function('y')

# The unevaluated Laplace transform:
Y = laplace_transform(y(t), t, p)

eq = Eq(diff(y(t), (t, 2)) + 2 * diff(y(t), t) + 16*y(t), 
        cos(4*t))

The right-hand side looks like this:

But how do we describe

in sympy ? Use the Subs class! It represents unevaluated substitutions of an expression. That's exactly what we need here. So our initial conditions are

Now we can write the Laplace transform of the differential equation as

Solve this for L[y]:

and transform back from Laplace to normal space:

Tidy up a little bit:

And this is a neat form where we will stop. The Heaviside 𝜃θ-function makes all values with 𝑡<0 equal to zero, which is ok, since we only need solutions for 𝑡≥0.

Note that the Laplace method automatically takes care of the initial conditions without the pain of determining constants from a general solution! That makes it a lot more comfortable than most other methods.

Finally, let’s plot the solution:

As you can see, the solution (blue) needs some time until it follows the external force. Finally, it oscillates with the same frequency, but with a constant phase shift.

The example above showed how to solve the problem with homogeneous initial conditions (𝑦(0)=𝑦′(0)=0). But using the Laplace technique is certainly not restricted to that. Simply plug in the nonhomogeneous initial conditions, solve for 𝑌(𝑝), do the inverse Laplace transform, and, voilà, you’ve got the solution 𝑦(𝑡).