I want to solve a coupled system of ODEs in matrix form 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
in matrix form and 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
time=np.linspace(0,4,1/dt)
y_init=np.ones(2)*-15.
y_tot=odeint(dy_dt,[y_init,m_init],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_tot=odeint(dy_dt,[y_init,m_init],time)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/odepack.py", line 215, in odeint
ixpr, mxstep, mxhnil, mxordn, mxords)
ValueError: Initial condition y0 must be one-dimensional.
The initial value to odeint must be an array, not a matrix. Try use y0=np.hstack((y_init, m_init)) and put that as the initial value (y0 is the second argument to odeint).