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

我试图解决我的论文在化学微分方程，有我绊了关于SciPy的的微分方程求解“odeint”的问题。SciPy的：实施微分方程的两种方法：两种不同的解决方案：回答

首先我由函数CIDNP_1实现差分（CIDNP是一个化学现象，这解释了不寻常的变量）根据SciPy的网站上的例子。但是，即使正确的方向，解决方案也是如此。

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给出错误的结果任何人都可以给我一个提示，我会为那么至少最后两个小时内未完全丧失很幸运。

问候，

雅各布

来源

2017-09-06 JakobJakobson

您确认输入'y'是语义上的时间导数，而不仅仅是国家回答'y = [P，Q]'RESP。 'P，Q = y'？这种重新贴标签也会避免发生混淆。 – LutzL

Semantivally你是对的。微分方程描述了系统的状态并且没有导数。但是，因为我最初只是想在五分钟内编码，所以我没有太多考虑。在下一个版本中，我会改变它。 – JakobJakobson

在你的第一个版本，你执行更新不是在同一时间，当你执行到线

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

不simulatnously;因此，您使用已经更新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]

这会解决您的问题。

来源

2017-09-06 20:06:14 FlashTek

非常感谢。下次我会尽量不在同一行中覆盖我的变量。 Regards – JakobJakobson

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

来源

2017-09-06 20:12:19

非常感谢。下次我会尽量不在同一行中覆盖我的变量。问候 – JakobJakobson

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

回答

相关问题