Python has many great inbuilt packages that make the 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 equationsa,b,c = sp.var(‘a b c’)#Now define expressions on these variablesExpression = 2*a*b + b*c – 8*a*b*c#We can find the variables used in an expression by .atoms() functionSymbols = Expression.atoms(sp.Symbol)Constants = Expression.atoms(sp.Integer)#We can substitute the values of variables used in an expression partially or completelyMod_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 functionsFunc = sp.lambdify(Symbols,Expression,”numpy”)#Once converted they can be used just like normal lambda functionsRes = Func(1,2,3)[/code]

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 spimport numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import ode#Define variables and expressionsy1, y2, y3 = sp.var(‘y1 y2 y3’)exp1 = -0.04*y1 + 10*y2*y3exp2 = 0.04*y1- 10**4*y2*y3 – 3*10**7*y2**2exp3 = 3*10**7*y2**2symbols = [y1,y2,y3]#Convert the expressions into lambda functionsfunc1 = sp.lambdify(symbols,exp1,”numpy”)func2 = sp.lambdify(symbols,exp2,”numpy”)func3 = sp.lambdify(symbols,exp3,”numpy”)fun=[func1,func2,func3]#Now define the step function(input to the scipy.integrate.ode function)def int_step(t,y):result=[]for i in range(len(y)):result.append(fun[i](*y))return result#Here we define the integratordef int_ode(y0,time):res=[]res.append(y0)r = ode(int_step)r.set_integrator(‘vode’, method=’bdf’)r.set_initial_value(y0,time[0])for t in time[1:]:r.integrate(t)res.append(r.y)return res#Define time list and  initial valuetime = [10**i for i in range(4)] +[0] + [-10**j for j in range(4)]time = sorted(time)y0=[1,0,0]#Calculate the resultr = int_ode(y0,time)#plot y2 with respect to timeplot_data= [k[1] for k in r]plt.plot(time, plot_data)plt.savefig(‘y2.svg’)[/code]

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

For a 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 spa,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[/code]

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.

Check the amazing data analysis capability of Polly, book a demo today!