如何在python中用2个变量（+时间）求解2个微分方程

Question

im在python中非常新，我一直在尝试求解具有2个未知数的2个联立微分方程组。两个变量T,X都是时间的函数，T的导数：dT/dt取决于T和X，而dX_dt仅取决于X。我的第一个想法是创建一个函数fun(t,T,X)，该函数返回的值dT/dt和dX/dt并使用solve_ivp，但效果不佳。有人可以告诉我解决这类ODE系统的正确方法是什么？提前致谢。这是脚本：

# Imports
import math
import numpy as np
from scipy.integrate import solve_ivp

# PARTICLE 

R =50*10**(-6) # Particle Radius (m)
rho = 575.4 # Initial Particle density  (kg/m^3)
mp0 = rho*4*math.pi*R**3/3 # Particle's initial mass
#t_c = 0.60 # C mass fraction in coal
t_ash = 0.245 #initial ash mass fraction in coal
mc0=mp0*(1-t_ash) # initial carbon mass fraction in coal
SSA =100000 # Coal Specific Surface Area (kg/m^2)
ep = 0.85 # Particle's emissivity
E= 153299 # Carbon combustion activation energy (J/mole)

# WALL AND AIR TEMPERATURES and pressure

Tw = 1000 # Wall temperature (K)
Ta = 300  # Air temperature (K)
p = 101325 # Aur pressure (Pa)


# OTHER CONSTANTS

sb = 5.670374419*10**(-8) # Stephan-Boltzman Constant ()
Rg = 8.31446261815324 # Gas constant (J/K mol)
Mc = 12.0107*10**(-3) # Carbon Molar Mass (kg/mole)
X_O2 = 0.20947 # Oxygen's Volumetric fraction in Air
Co2inf = X_O2*p/(Rg*Ta) # # Oxygen molar concentration in air (moles/m^3)

def fun(t,T, X):

    # Cp of carbon (J/kg K)
    Cpc = 3*Rg/Mc*math.exp(1200/T)*((math.exp(1200/T)-1)/(1200/T))**2
    # Cp of ash (J/kg K)
    Cpash = 754+0.586*(T-273.15)
    # Reaction Enthalpy of CO (J/mole -absolute value - temperature correction)
    DHco = 110.5*10**3 + 6.38*(T-298.15)
    # Reaction Enthalpy of CO2 (J/mole -absolute value - temperature correction)
    DHco2 = 393.5*10**3 + 13.1*(T-298.15)
    # CO/CO2 ratio
    CO_CO2 = 10**2.5*math.exp(-25000/(Rg*T))
    # Particle-Air boundary layer mean temperature
    Tm = (T+Ta)/2
    # Air thermal conductivity (W/m K -was used for gas with different O2 content!!!!)
    lamda = 0.0207*(Tm/273.15)**2*(1+113.15/273.15)/(1+113.15/Tm)
    # Convection coefficient (W/m^2 K - Nusselt=2)
    h = lamda/R
    # Effective diffusivity of oxygen within the particle, (m^2/s)
    Do2 = 1.78*10**(-5)*(Tm/273.15)**1.75
    # Mass transfer coefficient of oxygen (m/s - Sherwood=2)
    kd = Do2/R
    # Kinetic constant (m/s)
    ks = 1050*math.exp(-153200/(Rg*T))
    # L parameter
    L = (CO_CO2+1)/(CO_CO2/2+1)
    # Oxygen Flux
    FluxO2 = Co2inf/(3*L/(ks*rho*4*math.pi*R**3*SSA)+1/(kd*4*math.pi*R**2))

    # HEAT BALANCE EQUATION
    dT_dt= (4*math.pi*R**2*ep*sb*(Tw**4-T**4) +  4*math.pi*R**2*h*(Ta-T) + (CO_CO2*DHco+DHco2)*FluxO2/(1+CO_CO2/2))/(mc0*(1-X)*Cpc+mp0*t_ash*Cpash)

    # Combustion Reaction
    dX_dt= Co2inf*Mc*ks*SSA*(1-X)

    return [dT_dt, dX_dt]

Init_cond=[300,0]

OutputTimes = np.linspace(0, 20, 100)
ans = solve_ivp(fun, (0, 20), Init_cond, method='RK45', t_eval=OutputTimes)

Answer 1

返回值（状态的导数）是元组（或概念上的向量）的方式相同，ODE函数还必须接受状态作为相同维的元组/向量。该求解器基于通过向量值向量函数实现的ODE系统。

def fun(t,u):
    T, X = u;
    # etc.

通过此更改，您的程序应运行或至少给出更多有趣的错误。

---

[通常情况下，您可能想尝试method="Radau"是否能以更少的内部步骤实现更快的集成。