通过 Python 中的矩阵形式求解两组耦合 ODE

Question

我想求解具有以下形式的两组变量（即 {y} 和 {m}）的矩阵形式的 ODE 耦合系统：

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

其中 C 是一个矩阵，[2 1, -1 3]。

另一方面，我想解这些方程：

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

最终能够通过矩阵形式绘制 ys 和 ms 与时间的关系图。我写了下面的程序：

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

但是我遇到了以下错误：

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

谁能帮帮我！

Answer 1

看看这个以了解发生了什么：

# we have a collection of different containers, namely list, tuple, set & dictionary
master = [[1, 2], (1, 2), {1, 2}, {1: 'a', 2: 'b'}]

for container in master:
    a, b = container  # python will automatically try to unpack the container to supply a & b with values
    print(a, b)  # all print: 1 2 since a = 1 and b = 2 after the unpacking

如果我有一个容器，其值多于我尝试提供的变量，我会收到 "too many values to unpack" 错误，例如：

container = [1, 2, 3]
a, b = container  # this raises an error, the value 3 has nowhere to go

但是您可以通过以下方式说 "dump all the rest to b"：

a, *b = container 
print(a, b)  # -> 1 [2, 3] so a = 1 and b = [2, 3]

回到你知道的，当你说：y, m = Y 你必须确保 Y 是一个恰好有 2 个对象的容器，但事实并非如此。最后，正如我在评论中所说，您似乎没有在任何地方调用函数 dy_dt。

Answer 2

odeint 适用于一维数组。在你的函数 dy_dt 中，Y 将作为长度为 4 的一维 numpy 数组传入。要将其拆分为 y 和 m，你可以编写

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

函数 dy_dt 还必须 return 一个长度为 4 的序列，包含 [y[0]', y[1]', m[0]', m[1]' ]（即四个量的导数）。看起来你的函数目前 return 在 dy 中只有两个量，所以你也必须修复它。

这一行

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

看起来很像欧拉的方法。这不是你想要的。根据您在问题中所写的内容，您可以这样写：

dm = -4. * m * y

那么你可以 return np.hstack((dy, dm)).

总结一下，你可以这样写dy_dt：

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]]

通过 Python 中的矩阵形式求解两组耦合 ODE

Solving two sets of coupled ODEs via matrix form in Python

python

numpy

scipy

odeint