如何使用odeint来差分ode解决方案

如何解决如何使用odeint来差分ode解决方案

我正在使用odeint来解决python中的微分方程，我想计算ode解的导数。在MATLAB中，ode解的导数可以计算如下：

ode_fcn = @(t,y) ...;                                       % ODE Function
[t,y] = ode15i(ode_fcn,...);                               % Integrate ODE
dy = ode_fcn(t,y);                                          % Calculate Derivatives

但是，在python中，我尝试使用以这种方式定义的ode函数，这会导致错误的答案。例如，我定义了dmdT函数，如下所示：

dmdT(m,T,*parameters):
    pass
return [dm1dT,dm2dT,...]

然后我解决了ode函数

ode_solution = odeint(dmdT,m0,T_span,args=parameters)

我希望按以下方式计算ode_solution的导数：

derivative = dmdT(ode_solution,*parameters)

但这有时可能是错误的，我不知道为什么。

所以我的问题是：我可以像MATLAB一样使用定义的ode函数来计算ode解的导数吗？

另一种方法是在Numpy中使用diff函数，但它将数据长度缩短1。因此，我使用中差法定义了一个新的diff函数，该函数不同于Numpy中的diff来解决此问题。但是这种方法还不能像MATLAB中那样优雅。

我定义的diff函数如下：

def diff(x,y):
    x = np.array(x)
    y = np.array(y)
    # Forward difference
    dydx_start = (y[1]-y[0])/(x[1]-x[0])
    # Central difference quotient
    dydx = np.diff(y)/np.diff(x)
    dydx_mid = (dydx[1:] + dydx[:-1])/2
    # Backward difference
    dydx_end = (y[-1]-y[-2])/(x[-1]-x[-2])    
    dydx_mid = np.append(dydx_start,dydx_mid)
    return np.append(dydx_mid,dydx_end)

任何建议都可能有帮助，非常感谢！

---

根据Warren Weckesser的建议，我试图说明我提到的问题。这是我使用的代码

def dmdT(m,beta,A1,E1,n1,k1,A2,E2,n2,k2):
    R = 8.314
    """
    FR -> char + volitale
    polyester -> char + volitale
    """
    if len(m) == 4:
        m_FR= m[0]
        m_polyester = m[1]
        m_residue = m[2]
        m_total = m[3]    

    
        rxn1 = -m[0]*A1 * np.exp(-E1/R/T) * (m_FR/m[0])**n1 / beta
        rxn2 = -m[1]*A2 * np.exp(-E2/R/T) * (m_polyester/m[1])**n2 / beta
    else:
        m_FR= m[:,0]
        m_polyester = m[:,1]
        m_residue = m[:,2]  
        m_total = m[:,3]
        
        rxn1 = -m[0,0] * A1 * np.exp(-E1/R/T) * (m_FR/m[0,0])**n1 / beta 
        rxn2 = -m[0,1] * A2 * np.exp(-E2/R/T) * (m_polyester/m[0,1])**n2 / beta 
    
    r1 = rxn1
    r2 = rxn2
    r3 = -k1 * rxn1 - k2 * rxn2

    
    dm_FRdT = r1
    dm_polyesterdT = r2
    dm_residuedT = r3
    dm_totalT = r1 + r2 + r3 
    
    return [dm_FRdT,dm_polyesterdT,dm_residuedT,dm_totalT]    

m0 = [0.0954577747,0.9045422253,1]
T = np.arange(300,1000,0.1)
parameters = [5,9.42351217e+04,5.08266345e+04,9.87252218e-01,6.90982451e-01,8.28675597e+16,2.29312948e+05,6.94773512e-01,4.00782605e-01]
ode_solution = odeint(dmdT,args = tuple(parameters))

MLR = dmdT(ode_solution,*parameters)[-1]
MLR2 = diff(T,ode_solution[:,-1])
plt.figure()
plt.plot(T,MLR)
plt.plot(T,MLR2)

这两种方法的解决方案是不同的。当我使用dmdT求导数时，出现警告：

__main__:22: RuntimeWarning: invalid value encountered in power
__main__:23: RuntimeWarning: invalid value encountered in power

所以我想知道为什么会发生以及如何解决此问题？ odeint是否给出错误答案？

谢谢！