我正在尝试解决涉及使用常微分方程组的固定床反应器的问题，但由于某种原因它陷入了无限循环。这是代码：
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
" Dependant variables : Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P"
def f(W, y):
Fet, Fac, Fo, Fp, Fa, Fc, Ft, Pet, Pac, Po, Pp, Pa, Pc, P = y
r1 = 0.1036*np.exp(-3674/T)*((Po*Pet*Pac*(1+1.7*Pa)))/((1+0.583*Po*(1+1.7*Pa))*(1+6.8*Pac))
r2 = 1.9365*(10**5)*np.exp(-10116/T)*(Po*(1+0.68*Pa))/(1+0.76*Po*(1+0.68*Pa))
m = Fet * 28.06 + Fac * 60.052 + Fo * 32 + Fp * 86.09 + Fa * 18.02 + Fc * 44.01 # [ lb / min ]
G = m / A_c # [ lb / (min * ft ^2) ]
beta = (1.75 * G**2 / (rho_0 * g_c * D_p)) * ((1-phi) / phi**3)
dFet_dW = -r1 - r2
dFac_dW = -r1
dFo_dW = -r1/2 - 3*r2
dFp_dW = r1
dFa_dW = r1 + 2*r2
dFc_dW = 2*r2
dFt_dW = -r1/2
dP_dW = -(beta/(A_c*(1-phi)*rho_c))*(P_0/P)*(Ft/Ft_0)
dPet_dW = (Ft*P*dFet_dW-Fet*P*dFt_dW+Fet*Ft*dP_dW)/Ft**2
dPac_dW = (Ft*P*dFac_dW-Fac*P*dFt_dW+Fac*Ft*dP_dW)/Ft**2
dPo_dW = (Ft*P*dFo_dW -Fo*P*dFt_dW +Fo*Ft*dP_dW)/Ft**2
dPp_dW = (Fp*P*dFp_dW -Fp*P*dFt_dW +Fp*Ft*dP_dW)/Ft**2
dPa_dW = (Fa*P*dFa_dW -Fa*P*dFt_dW +Fa*Ft*dP_dW)/Ft**2
dPc_dW = (Fc*P*dFc_dW -Fc*P*dFc_dW +Fc*Ft*dP_dW)/Ft**2
return np.array([dFet_dW, dFac_dW, dFo_dW, dFp_dW, dFa_dW, dFc_dW, dFt_dW, dPet_dW, dPac_dW, dPo_dW, dPp_dW, dPa_dW, dPc_dW, dP_dW])
# Integration points
W_span = np.array([0, 100])
points = np.linspace(W_span[0], W_span[1], num = 101)
# Initial conditions
P_0 = 128 # psi
T = 423 # Kelvin
Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0 = 11, 11, 1.9, 0, 0, 0
Ft_0 = Fet_0 + Fac_0 + Fo_0
Pp_0, Pa_0, Pc_0 = 0, 0, 0
Pet_0 = (Fet_0/Ft_0)*P_0
Pac_0 = (Fac_0/Ft_0)*P_0
Po_0 = (Fo_0/Ft_0)*P_0
P_0 = Pet_0 + Pac_0 + Po_0
y_0 = np.array([Fet_0, Fac_0, Fo_0, Fp_0, Fa_0, Fc_0, Ft_0, Pet_0, Pac_0, Po_0, Pp_0, Pa_0, Pc_0, P_0])
# Constants
phi = 0.5 # Void fraction
A_c = 0.0122718463 # straight section area of the catalyst [ ft ^2 ]
rho_c = 60 # dry catalyst density [ lb / ft ^3 ]
g_c = 115826.4 # gravity conversion factor [ lb * ft / (min ^2 *lbf) ]
D_p = 0.020833333 # average particle diameter [ ft ]
M_0 = (28.05 * Fet_0 + 60.052 * Fac_0 + 31.998 * Fo_0) / (Ft_0) # molar mass of the initial mixture [ lb / lbmol]
rho_0 = (128 * M_0) / (10.73 + 731.07) # Initial density [ lb / ft ^3 ]
# Solution IVP
sol = solve_ivp(f, W_span, y_0, t_eval=points)
W = sol.t
Fet = sol.y[0]
Fac = sol.y[1]
Fo = sol.y[2]
Fp = sol.y[3]
Fa = sol.y[4]
Fc = sol.y[5]
Ft = sol.y[6]
Pet = sol.y[7]
Pac = sol.y[8]
Po = sol.y[9]
Pp = sol.y[10]
Pa = sol.y[11]
Pc = sol.y[12]
P = sol.y[13]
我尝试在代码中编写一些 print() ，同时试图了解可能发生的情况，并意识到一些值，如 Fet、Fac、Fo、Fp、Fa、Fc、Ft、m、G 和测试版没有改变。另一方面，r1、r2、Pet、Pac、Po、Pp、Pa、Pc和P周期性变化。
的代码似乎进入了无限循环，并非因为
solve_ivp
本身存在问题，而是因为的微分方程的定义方式导致了数值积分中的振荡或僵硬行为。具体来说，问题出在计算偏导数的方式上，例如
dPet_dW
、
dPac_dW
等。这些方程中包含了乘法和除法的组合，如果分母很小或接近于零，则会导致数值不稳定和振荡。
为了解决这个问题，可以尝试以下几种方法：
-
重新排列微分方程： 尽可能地重新排列的微分方程，以避免除法或使用分母很小的项。例如，对于
dPet_dW，可以尝试对其进行代数化简，使其不包含Ft作为分母。可以尝试对其他偏导数方程执行类似的操作。 -
使用隐式求解器：
solve_ivp默认使用显式求解器，对于僵硬方程组来说可能不稳定。尝试使用隐式求解器，例如 'Radau' 或 'BDF'。可以通过将method参数传递给solve_ivp来实现，例如：python sol = solve_ivp(f, W_span, y_0, t_eval=points, method='Radau') -
调整求解器公差：
solve_ivp使用绝对和相对公差来控制求解的精度。尝试减小这些公差以查看是否有帮助。可以通过将atol和rtol参数传递给solve_ivp来实现，例如：python sol = solve_ivp(f, W_span, y_0, t_eval=points, atol=1e-8, rtol=1e-6) -
检查初始条件和参数： 确保的初始条件和参数值合理且在预期范围内。错误的值可能会导致数值不稳定性。
-
简化模型： 如果可能，尝试简化的模型以减少方程的数量或复杂性。这可能有助于识别问题的根源并更容易找到解决方案。
总的来说，解决数值积分中的无限循环或僵硬行为需要仔细分析的方程和数值方法。尝试上述方法并观察它们如何影响的解决方案。如果问题仍然存在，请提供有关的方程、参数和初始条件的更多详细信息，以便我可以提供更具体的帮助。
标签：python,infinite-loop,ode From： 78776670