首页 文章

Scipy：实现微分方程的两种方法：两种不同的解决方案：回答

提问于
浏览 次
1

我试图解决我的化学论文的微分方程，在那里我偶然发现了一个关于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会给出错误的结果，我会非常幸运，因为至少在最后两个小时内并没有完全丢失 .

问候，

雅各布

python numpy scipy differential-equations

2 回答

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

    这可以解决你的问题 .

    回复于
  • 1

    在 CIDNP_1 中，在使用新值计算 dQ_dt 之前更改 dP_dt 的值：

    dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2   # changed dP_dt!
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2   # use the changed value!

    在 CIDNP_2 中，您可以同时计算它们，即 dQ_dt 使用 dP_dt 的原始值计算，而不是更改的值 . 你可以这样想

    a = -kt*dP_dt*R(t) - kt*beta*(R(t))**2       # no overwriting -
b = +kt*dP_dt*R(t) + kt*beta*(R(t))**2       # uses original value of dP_dt
return [a, b]

    你也可以加快速度

    def CIDNP_3(y, t):
    dP_dt, dQ_dt = y
    R_t = R0 / (1 + kt * R0 * t)
    res = kt * R_t * (dP_dt + beta * R_t)
    return [-res, res]
    回复于

相关问题