Solving Second Order Differential Equations in Python

This presentation outlines how to solve second order differential equations in python. The solution is obtained numerically using the python SciPy ode engine (integrate module). This solution is therefore not in analytic form but the output is as if the analytic function was computed for each time step.

The method used here is generally applicable to solving higher order differential equations in python as well.

We will use a series RLC circuit as our second order differential equations example.

Step 1

Split your second order differential equation into two first order ordinary differential equations (ODE).

Example, if the second order differential equation is:

            ax" + bx' + c = k                (1)

we can let:

            q = x'                    (2)

if we substitute this value in the original equation we now have:

            aq' + bq + c = k              (3)

Equations (2) and (3) form a 'system' of first order ordinary differential equations that represent the second order equation shown in equation (1).

For our series RLC circuit we have:

            LCv" + RCv' + v = Vs            (4)
where v is the voltage across the capacitor

Equation (4) is a second order differential equation in the capacitor voltage.

If we let:

            x = v'                  (5)

we then substitute this variable into equation (4) to obtain:

            LCx' + RCx + v = Vs              (6)

We now have two first order differential equations to represent our circuit.

Step 2

Now we want to represent out system of differential equations in python. First we rearrange the first order differential equations in the form:

            x' = f(x)                (7)

Applying this to our series RLC circuit example we get:

            v' = x                  (8)
            x' = (Vs - v - RCx)/LC           (9)

Step 3

Now we can write some code.

Implement a python function that returns the right hand sides of the two first order differential equations.

For our example the function is:

def rlc(A,t): 
    Vc,x=A
    V = 1.0 #voltageSource
    R = 5.0
    L=100.0e-9 #100nH
    C = 1.0e-9 #1nF
    res=array([x,(V-Vc-(x*R*C))/(L*C)])
    return res
            

Now to solve this we have to define our time step and duration. We can do this using the linspace function as follows:

time = linspace(0.0,0.6e-6,1001)

We cane now call the scipy integrate function to obtain the solution and transpose it to get two columns with the values at each time step. This is done as follows:

vc,x = integrate.odeint(rlc,[0,0],time).T

Step 4

We can now plot our solution for easy visualization:

figure()
plot(time,vc)
xlabel('t')
ylabel('Vc')
show()
        

This gives:

That's it. Don't forget to share this post with anyone you think might find it useful.

Happy coding 🙂

Full Code

# -*- coding: utf-8 -*-
from scipy import integrate
from pylab import * # for plotting commands

def rlc(A,t):
    Vc,x=A
    V = 1.0 #voltageSource
    R = 5.0
    L=100.0e-9 #100nH
    C = 1.0e-9 #1nF
    res=array([x,(V-Vc-(x*R*C))/(L*C)])
    return res

time = linspace(0.0,0.6e-6,1001)
vc,x = integrate.odeint(rlc,[0.0,0.0],time).T
figure()
plot(time,vc)
xlabel('t')
ylabel('Vc')
show()