无法使用solve_ivp得到正确的解决方案

Question

我正在尝试使用 solve_ivp 求解微分方程，但我没有得到正确的解决方案。但是，我使用 ideint 获得了正确的解决方案。我的solve_ipv程序有问题吗？

ODEINT 程序：

# Arhenius Function
def Arhenius(a, T):
    dadT = np.exp(lnA)/v * np.exp(- E / (8.3144 * T)) * c * np.abs(a) ** m * np.abs((1 - np.abs(a))) ** n
    return dadT

# Initial data
pt = 100000
T0 = 273
Tf = 1500
a0 = 0.0000000001
T = np.linspace(T0, Tf, pt)

a_sol = np.zeros(pt)
dadt = np.zeros(pt)

# ODE solve
a_t = odeint(Arhenius, a0, T)

# For removing errored values and have maximum at 1 
search1 = np.where(np.isclose(a_t, 1))
try:
    ia_1 = search1[0][0]
    a_sol[0,:] = a_t[:,0]
    a_sol[0,ia_1+1:pt] = 1
except:
    a_sol = a_t[:,0]

#  Calculate the new derivative
dadt = np.exp(lnA) * np.exp(- E / (8.3144 * T)) * c * np.abs(a_sol) ** m * np.abs((1 - a_sol)) ** n

使用solve_ivp编程：

# Arhenius Function
def Arhenius(a, T):
    dadT = np.exp(lnA) * np.exp(- E / (8.3144 * T)) * c * np.abs(a) ** m * np.abs((1 - np.abs(a))) ** n
    return dadT

# Initial data
pt = 100000
T0 = 273
Tf = 1500
a0 = 0.0000000001
T = np.linspace(T0, Tf, pt)

a_sol = np.zeros(pt)
dadt = np.zeros(pt)

# ODE solve
a_t = solve_ivp(Arhenius, t_span = (T0, Tf), y0 = (a0,), t_eval = T, method = 'RK45')
a_sol= a_t.y

#  Calculate the new derivative    
dadt = np.exp(lnA) * np.exp(- E / (8.3144 * T)) * c * np.abs(a_sol) ** m * np.abs((1 - a_sol)) ** n

Answer 1

solve_ivp 将 a 和 T 颠倒。还需要进行一些修复。代码现在的工作原理如下：

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

lnA = 36
E = 125000
c= 1
m = 0
n= 1.5
v = 5

# Arhenius Function
def Arhenius(T, a):
    if type(a) is not float:
        a = a[0]
    dadT = np.exp(lnA) * np.exp(- E / (8.3144 * T)) * c * np.abs(a) ** m * np.abs((1 - np.abs(a))) ** n
    return [dadT]

# Initial data
pt = 100000
T0 = 273
Tf = 450
a0 = 0.0000000001
T = np.linspace(T0, Tf, pt)

# ODE solve
a_t = solve_ivp(Arhenius, t_span = (T0, Tf), y0 = (a0,), t_eval = T, method = 'Radau')
a_sol = a_t.y[0]
T = a_t.t


#  Calculate the new derivative    
dadt = np.exp(lnA) * np.exp(- E / (8.3144 * T)) * c * np.abs(a_sol) ** m * np.abs((1 - a_sol)) ** n

plt.plot(a_t.t, a_t.y[0])
plt.xlim([0, Tf])