我的odeint实现是否存在逻辑错误？

Question

我得到了错误的解决方案，不确定odeint是否是解决此ODE系统的正确工具。

我正在尝试通过求解ODE系统来对简单的一阶化学反应进行建模。从逻辑的角度来看，我的函数是正确的，我可以在MATLAB中解决这个问题，几乎没有问题。我也希望能够在python中完成这项工作。我认为odeint是完成这项工作的工具，但我可能是错的。我的解决方案不应每次都收敛于自变量= 10，但无论输入如何，它始终会收敛。

from matplotlib.pyplot import (plot,grid,xlabel,ylabel,show,legend)
import numpy as np 
from scipy.integrate import odeint

wght= np.linspace(0,20)
# reaction is A -> B
def PBR(fun,W):

    X,y = fun

    P_0=20;#%bar
    v_0=5; #%m^3/min
    y_A0=1; #unitless
    k=.005; #m^3/kg/min
    alpha=0.1; #1/kg
    epi=.13; #unitless
    R=8.314; #J/mol/K

    F_A0= .5 ;#mol/min

    ra = -k *y*(1-X)/(1+epi*X)


    dX = (-ra)/F_A0
    dy = -alpha*(1+epi*X)/(2*y)

    return [dX,dy]

X0 = 0.0
y0 = 1.0


sol = odeint(PBR, [X0, y0],wght)


plot(wght, sol[:, 0], 'b', label='X')
plot(wght, sol[:, 1], 'g', label='y')
legend(loc='best')
xlabel('W')
grid()
show()
print(sol)

输出图

Answer 1

您的图形不可重现，在python 3.7，scipy 1.4.1中，我在12.5左右达到几乎无穷大，并且在清理工作空间后，两个图形在时间10之后都移为零。

问题是您在第二个等式中除以2*y。实际上，这意味着y是某些其他函数的平方根，并且该函数的类型在正切值变为垂直的情况下在较小的数值上表现不佳，此后不久便变为零。

但是，可以用一种简单的方法将除法分解为整数

    dy = -alpha*(1+epi*X)*y/(1e-10+2*y*y)

或者以更复杂的方式使解决方案远离零，实际上并没有改变

    dy = -alpha*(1+epi*X)*y/(y*y+max(1e-12,y*y))

结果图

如果仍然不是预期的行为，则相对于Matlab版本，您的ODE系统中会有其他差异。