Solving two sets of coupled ODEs via matrix form in Python

Question

I want to solve a coupled system of ODEs in matrix form for two sets of variable (i.e. {y} and {m}) which has such a form:

y'_n = ((m_n)**2) * y_n+(C * y)_n , m'_n=-4*m_n*y_n

where C is a matrix, [2 1, -1 3].

On the other hand I want to solve these equations:

y'1= m1 ** 2 * y1 + 2 * y1 + y2
y'2= m2 ** 2 * y2 - y1 + 3 * y3
m'1= -4 * m1 * y1 ,
m'2= -4 * m2 * y2
y1(0)=y2(0)=-15. and m1(0)=m2(0)=0.01

to finally be able to plot ys and ms versus time via matrix form. I wrote the following program:

import numpy as np
from pylab import plot,show
from scipy.integrate import odeint

C=np.array([[2,1],[-1,3]])
dt=0.001
def dy_dt(Y,time):
 y,m=Y
 m=m+dt*(-4.*m*y)
 dy=m**2*y+np.dot(C,y)
 return dy

m_init=np.ones(2)*0.01
y_init=np.ones(2)*-15.
time=np.linspace(0,4,1/dt)
y0=np.hstack((y_init, m_init))
y_tot=odeint(dy_dt,y0,time)



plot(time,y_tot[0])#y_1
plot(time,y_tot[1])#y_2
plot(time,y_tot[2])#m_1
plot(time,y_tot[3])#m_2
show()

but I encountered the following error:

 y,m=Y   
 ValueError: too many values to unpack 

Can anybody help me!

Answer 1

odeint works with one-dimensional arrays. In your function dy_dt, Y will be passed in as a one-dimensional numpy array with length 4. To split it into y and m, you can write

y = Y[:2]
m = Y[2:]

The function dy_dt must also return a sequence of length 4, containing [y[0]', y[1]', m[0]', m[1]'] (i.e. the derivatives of the four quantities). It looks like your function currently returns just two quantities in dy, so you will have to fix that, too.

This line

m=m+dt*(-4.*m*y)

looks suspiciously like Euler's method. That's not what you want here. Based on what you wrote in the question, you could write:

dm = -4. * m * y

Then you could return np.hstack((dy, dm)).

To summarize, you can write dy_dt like this:

def dy_dt(Y, time):
    y = Y[:2]
    m = Y[2:]
    dy = m**2 * y + np.dot(C, y)
    dm = -4 * m * y
    return np.hstack((dy, dm))  # or:  return [dy[0], dy[1], dm[0], dm[1]]