python - Scipy:实现微分方程的两种方法:两种不同的解决方案:已回答

Question

我试图为我的化学论文求解一个微分方程，在那里我偶然发现了一个关于 scipy 的微分方程求解器“odeint”的问题。

首先，我根据scipy网站上的示例，通过函数CIDNP_1(CIDNP是一种化学现象，解释了不寻常的变量)实现了微分。但解决方案甚至偏离了正确的方向。

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

R0 = 5e+5
kt = 5e5/R0
beta = 3/R0

def CIDNP_1(y, t):
    dP_dt, dQ_dt = y

    def R(t):
        return R0/(1 + kt*R0*t)

    dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
    dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2
    return [dP_dt, dQ_dt]


def CIDNP_2(y, t):    
    dP_dt, dQ_dt = y

    def R(t):
        return R0/(1 + kt*R0*t)

    return [-kt*dP_dt*R(t) - kt*beta*(R(t))**2, \
            +kt*dP_dt*R(t) + kt*beta*(R(t))**2]

y0 = [-1, +1]
t = np.linspace(1e-9, 100e-6, 1e3)
sol_1 = scipy.integrate.odeint(CIDNP_1, y0, t)
sol_2 = scipy.integrate.odeint(CIDNP_2, y0, t)

然后我将解决方案更改为 CIDNP_2，得到了正确的结果，但在我看来，实现没有任何区别，因为变量 dP_dt 和 dQ_dt 在实现 CIDNP_1 中没有更改。

因此，任何人都可以给我一个提示，为什么 CIDNP_1 的实现会给出错误的结果，我会非常幸运，因为至少最后两个小时没有完全丢失。

问候，

雅各布

Answer 1

在您的第一个版本中，您在执行 to 行时不同时执行更新

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2

不同时；因此，您可以使用已更新的 dP_dt 来更新 dQ_dt。这是 ODE 系统的错误实现。你的第二种方法更好，因为它没有这种错误。您要么必须直接返回更新后的值，要么必须在计算新的 dQ_dt 之前将新计算的 dP_dt 值保存在另一个变量中。

new_dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2
new_dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2

return [new_dP_dt, new_dQ_dt]

这可以解决您的问题。