Any way to solve a system of coupled differential equations in python?

Question

I've been working with sympy and scipy, but can't find or figure out how to solve a system of coupled differential equations (non-linear, first-order).

So is there any way to solve coupled differential equations?

The equations are of the form:

V11'(s) = -12*v12(s)**2
v22'(s) = 12*v12(s)**2
v12'(s) = 6*v11(s)*v12(s) - 6*v12(s)*v22(s) - 36*v12(s)

with initial conditions for v11(s), v22(s), v12(s).

Orthman

Take a look at sage. It offers mathmatica-like functionality with python syntax. It might be able to solve diff eqs. – Tonsillectomy

Are you looking for an analytical solution, or a numerical solution? (You mentioned using sympy, so you might be hoping for an analytical solution, if there is one.) – Tiberius

@WarrenWeckesser A numerical solution, similar to the NDsolve for mathematica. – Orthman

It's first order, but this doesn't look like a linear system to me, as you have powers and products of the dependent variables. – Hawking

@Bitrex You're right, I mistakenly wrote linear rather than non-linear. Post has been updated. Good catch! – Orthman

I doubt this has a symbolic (analytic) solution, but even if it did, SymPy's dsolve doesn't yet support systems of ODEs. – Kraul

Answer 1

For the numerical solution of ODEs with scipy, see scipy.integrate.solve_ivp, scipy.integrate.odeint or scipy.integrate.ode.

Some examples are given in the SciPy Cookbook (scroll down to the section on "Ordinary Differential Equations").

Tiberius

The documentation seemed a bit difficult to understand, but the Cookbook you linked fitted just perfect for what I was trying to do, thank you! – Aldred

Answer 2

In addition to SciPy methods odeint and ode that were already mentioned, it now has solve_ivp which is newer and often more convenient. A complete example, encoding [v11, v22, v12] as an array v:

from scipy.integrate import solve_ivp
def rhs(s, v): 
    return [-12*v[2]**2, 12*v[2]**2, 6*v[0]*v[2] - 6*v[2]*v[1] - 36*v[2]]
res = solve_ivp(rhs, (0, 0.1), [2, 3, 4])

This solves the system on the interval (0, 0.1) with initial value [2, 3, 4]. The result has independent variable (s in your notation) as res.t:

array([ 0.        ,  0.01410735,  0.03114023,  0.04650042,  0.06204205,
        0.07758368,  0.0931253 ,  0.1       ])

These values were chosen automatically. One can provide t_eval to have the solution evaluated at desired points: for example, t_eval=np.linspace(0, 0.1).

The dependent variable (the function we are looking for) is in res.y:

array([[ 2.        ,  0.54560138,  0.2400736 ,  0.20555144,  0.2006393 ,
         0.19995753,  0.1998629 ,  0.1998538 ],
       [ 3.        ,  4.45439862,  4.7599264 ,  4.79444856,  4.7993607 ,
         4.80004247,  4.8001371 ,  4.8001462 ],
       [ 4.        ,  1.89500744,  0.65818761,  0.24868116,  0.09268216,
         0.0345318 ,  0.01286543,  0.00830872]])

With Matplotlib, this solution is plotted as plt.plot(res.t, res.y.T) (the plot would be smoother if I provided t_eval as mentioned).

Finally, if the system involved equations of order higher than 1, one would need to use reduction to a 1st order system.

Lysine