Using symbolic mathematics to solve System of ODEs in python

by Kailash Yadav (

Point of contact: Raaisa (

Python has many great inbuilt packages that make solving system of ODEs a piece of cake. In this post I will explain how we can use Sympy, Scipy, Numpy and some other libraries to solve a system of ODEs.

We will use the Robertson stiff system of odes in this blog-

$latex \large X’ = -0.04*X + 10^7*Y*Z$
$latex \large Y’ = 0.04*X – 10^7*Y*Z -3.10^7*Y^2$
$latex \large Z’ = 3*10^7*Y^2$

Sympy stands for symbolic mathematics library in python. It allows you to create and manipulate mathematical expressions easily.

Let’s have a look on the below python code snippet –

[code lang=python]
import sympy as sp

#define the variables to be used in equations
a,b,c = sp.var('a b c')

#Now define expressions on these variables
Expression = 2*a*b + b*c – 8*a*b*c

#We can find the variables used in an expression by .atoms() function
Symbols = Expression.atoms(sp.Symbol)
Constants = Expression.atoms(sp.Integer)

#We can substitute the values of variables used in an expression partially or completely
Mod_exp1 = Expression.xreplace({a:10,b:20})
Mod_exp2= Expression.subs({a:10,b:20})
Mod_exp3 = Expression.evalf(subs={a:10,b:20})

#We can also convert these expressions into lambda functions
Func = sp.lambdify(Symbols,Expression,"numpy")

#Once converted they can be used just like normal lambda functions
Res = Func(1,2,3)


Here we converted a sympy expression into lambda function which will then run at numpy speed instead of slower sympy speed.

There are many other useful tools in sympy which I encourage you to try to firm your grip on sympy.

Now let us look at how to solve a system of ODEs in python with sympy –

Here we will take y = (y1,y2,y3) to be the vector (X’,Y’,Z’) defined at the very end of this blog.We will use scipy.integrate.ode with Vode integrator and BDF method.
This is a stiff system of odes. Vode integrator with BDF method works quite good in general for stiff odes.

Let us first define the system of odes-

[code lang=python]
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode

#Define variables and expressions
y1, y2, y3 = sp.var('y1 y2 y3')
exp1 = -0.04*y1 + 10*y2*y3
exp2 = 0.04*y1- 10**4*y2*y3 – 3*10**7*y2**2
exp3 = 3*10**7*y2**2
symbols = [y1,y2,y3]

#Convert the expressions into lambda functions
func1 = sp.lambdify(symbols,exp1,"numpy")
func2 = sp.lambdify(symbols,exp2,"numpy")
func3 = sp.lambdify(symbols,exp3,"numpy")

#Now define the step function(input to the scipy.integrate.ode function)
def int_step(t,y):
for i in range(len(y)):
return result

#Here we define the integrator
def int_ode(y0,time):
r = ode(int_step)
r.set_integrator('vode', method='bdf')
for t in time[1:]:
return res

#Define time list and  initial value
time = [10**i for i in range(4)] +[0] + [-10**j for j in range(4)]
time = sorted(time)

#Calculate the result
r = int_ode(y0,time)

#plot y2 with respect to time
plot_data= [k[1] for k in r]
plt.plot(time, plot_data)


Here we have used sympy, numpy and scipy to integrate a non-linear system of odes.

For large system of odes, defining all the equations by hand is not feasible and that’s exactly the scenario where sympy comes to the rescue. We can create expressions (and hence functions) from strings too with the help of sympify function. For example –

[code lang=python]
import sympy as sp
a,b,c = sp.var('a b c')
string = "a*b*c"
expr = sp.sympify(string)
fun = sp.lambdify([a,b,c],expr)
assert fun(1,2,3) == 6.0

In its raw form, Sympy is not the fastest symbolic library but its simplicity beats everything that i have tried in python. We can try many things to speed up sympy such as trying cython on sympy code, using lambdify numpy, theano and ufuncify instead of evalf and subs.







Vipul Jain Written by:

Be First to Comment

Leave a Reply

Your email address will not be published. Required fields are marked *