**by Kailash Yadav (kailash.yadav@elucidata.io)**

**Point of contact: Raaisa (raaisa@elucidata.io)**

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)

[/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 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")

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 integrator

def 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 value

time = [10**i for i in range(4)] +[0] + [-10**j for j in range(4)]

time = sorted(time)

y0=[1,0,0]

#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)

plt.savefig('y2.svg')

[/code]

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

[/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.