常微分方程指数误差

Question

我开始学习 Python（这是我第一次尝试 Python 脚本）。我正在尝试求解，然后绘制一系列微分方程，这些微分方程显示了以下方程 h2o+co2 o2+ch2o 中的化合物浓度之间的关系。

但是当我运行脚本时：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

y=[3.0,1.0,3.0,1.0]
g=15 
i=25 
tspan= np.linspace(0,5,100)

# h2o+co2 <-> o2+ch2o
def chemanalysis(y,t):

    x =i*math.sin(math.pi*t/12)**2

    dydt= -x*(y[0]*y[1]) + g*(y[2]*y[3]) #co2
    dydt= -x*(y[0]*y[1]) + g*(y[2]*y[3]) #h2o
    dydt= x*(y[0]*y[1])  - g*(y[2]*y[3]) #o2
    dydt= x*(y[0]*y[1])  - g*(y[2]*y[3]) #ch2o


    return dydt


y=odeint(chemanalysis,1,tspan)

print(y)

我收到以下错误：

Traceback (most recent call last):
File "first.py", line 25, in <module>
        y=odeint(chemanalysis,1,tspan)
        ixpr, mxstep, mxhnil, mxordn, mxords)
      File "first.py", line 16, in chemanalysis
        dydt= -x*(y[0]*y[1]) + g*(y[2]*y[3]) #co2
    IndexError: index 1 is out of bounds for axis 0 with size 1

我知道这可能是一个非常直接的问题，因为我缺乏一般的 python/编程经验。我真的很感激任何关于我哪里出错或我可以改进这段代码的帮助。我在其他线程上阅读了一些其他类似的问题/错误消息，但是我自己无法解决问题。我知道目前没有附加到这个脚本的图形组件，我想我会等到先解决这个问题。

问候

Answer 1

你想在这里做两件事：

将y=odeint(chemanalysis, 1, tspan) 更改为y2=odeint(chemanalysis, y, tspan)

要使odeint 工作，您需要返回与y 大小相同的列表。

def chemanalysis(y,t):
    x =i*math.sin(math.pi*t/12)**2

    dydt = []

    dydt.append(-x*(y[0]*y[1]) + g*(y[2]*y[3])) #co2
    dydt.append(-x*(y[0]*y[1]) + g*(y[2]*y[3])) #h2o
    dydt.append(x*(y[0]*y[1])  - g*(y[2]*y[3])) #o2
    dydt.append(x*(y[0]*y[1])  - g*(y[2]*y[3])) #ch2o


    return dydt

我们可以让这个看起来更像 Python，但为了简洁而省略。

---

关于为什么需要在chemanalysis 中使用append 的更多解释：

dydt= -x*(y[0]*y[1]) + g*(y[2]*y[3]) #co2
dydt= -x*(y[0]*y[1]) + g*(y[2]*y[3]) #h2o
dydt= x*(y[0]*y[1])  - g*(y[2]*y[3]) #o2
dydt= x*(y[0]*y[1])  - g*(y[2]*y[3]) #ch2o

return dydt

上面只会返回 dydt 作为最后一个分配的值 #ch2o，你会丢失其他数据并每次重写 dydt。而是将其设为一个列表，odeint 无论如何都要求您这样做。