Настройка ODEINT для нескольких ODE

Question

Я новичок в Python, и я пытаюсь перенести код MATLAB на Python. Здесь, пытаясь решить несколько ODE в Python с помощью решателя ODEINT scipy, я сталкиваюсь с проблемой в настройке начальных условий для передачи функции ODE. Ниже приведен код:

# Import all packages here

import numpy as np
import cantera as ct
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math as mt

# Define functions
def RMGCanODE(IC,t,gas,p0,patm,NM):
    Ru = ct.gas_constant/1000 # Universal gas constant [J/mol.K]
    Wk = gas.molecular_weights # Molecular weights of species [g/mol]
    X0, T0 = IC
    gas.TPX = T0,p0,X0 # Set the temperature, pressure, and mole fractions [K, atm, -]
    Mmix = gas.mean_molecular_weight # Mean molar mass of gas mix [kg/kmol]
    rhom = p0/Ru/T0/1e6 # Molar density [mol/cm^3]
    rhoM = p0*Mmix/Ru/T0/1000 # Mass density [kg/m^3]

    # Mass fraction calculations [-]
    Y0 = np.zeros([NM])
    for i in range(NM):
        Y0[i] = X0[i]*Wk[i]/Mmix

    # Thermodynamic property calculation [molar basis]
    Cpk = gas.standard_cp_R*Ru # Specific enthalpy at constant pressure for all species [J/mol.K]
    Hk = gas.standard_enthalpies_RT*Ru*T0 # Enthalpy of each species [J/mol]
    Sk = gas.standard_entropies_R*Ru # Entropy of each species [J/mol.K]

    # Mixture averaged specific heat [J/mol.K]
    Cpmix = 0
    for i in range(NM):
        Cpmix = Cpmix + Cpk[i]*X0[i]
    CpmixM = 1000*Cpmix/Mmix # Mixture averaged specific heat [J/kg.K]

    # Enthalpy of species at standard state [J/kg]
    HkM = np.zeros([NM])
    for i in range(NM):
        HkM[i] = Hk[i]/Wk[i]*1000

    # Mixture averaged entropy [J/mol.K]
    Skmix = 0
    for i in range(NM):
        if X0[i] != 0:
            Skmix = Skmix + X0[i]*(Sk[i] - Ru*(mt.log10(X0[i]*p0/patm))) # Entropy of species [J/mol.K]

    # Mixture averaged Gibbs free energy [J/mol]
    Gkmix = 0
    for i in range(NM):
        if X0[i] != 0:
            Gkmix = Gkmix + X0[i]*(Hk[i] - T0*(Sk[i] - Ru*(mt.log10(X0[i]*p0/patm))))

    # Call Cantera function for WdotM values [kg/m^3.s]
    Wdot = gas.net_production_rates # Net production rate of species [kmol/m^3.s]
    WdotM = np.zeros(NM)
    for i in range(NM):
        WdotM[i] = Wdot[i]*Wk[i]

    print(WdotM)

    # Mass fraction ODE
    dYdt = WdotM/rhoM

    # Energy ODE
    sum1 = 0
    for i in range(NM):
        sum1 = sum1 + HkM[i]*WdotM[i]/rhoM/CpmixM

    dTdt = sum1 # Energy equation, temperature change with time [K/s]

    # Mole fraction ODE
    dXdt = np.zeros(NM)
    for j in range(NM):
        sum2 = 0
        sum3 = 0
        for i in range(NM):
            sum2 = sum2 + Y0[i]/Wk[i]
            sum3 = sum3 + dYdt[i]/Wk[i]
        numer = dYdt[j]/Wk[j]*sum2 - Y0[j]/Wk[j]*sum3
        denom = sum2*sum2
        dXdt[j] = numer/denom

    dZdt = np.array([dXdt.tolist(),dTdt])
    return dZdt


# Driver
# Set all initial values here
p_exp = 3.02 # Experimental pressure value [atm]
patm = 101325 # Conversion: atm ---> Pa
p0 = p_exp*patm # [Pa]
T0 = 934 # Initial temperature [K]

# Set all switches here
writetofile = False # Toggle True or False to write to a .txt file
drawfig = True # Toggle True or False to plot data

# Solver parameters
tmin = 0 # Start time for simulation [s]
tmax = 3 # End time for simulation [s]
t0 = np.linspace(tmin,tmax,251) # Initial time conditions
atol = 1e-12 # Absolute tolerance
rtol = 1e-3 # Relative tolerance

# Reaction mechanism to use
filename = 'gri30.cti'
# Set up relevant data
gas = ct.Solution(filename) # Set the operating gas
gas.TPX = T0,p0,'H2:0.0095,O2:0.0049,N2:0.9856' # Set the temperature, pressure, and mole fractions
X0 = gas.X # Copy mole fractions in an array
Wk = gas.molecular_weights # Specify gas molecular weights
NM = gas.n_species # Number of species in the gas mix
v_p = ct.Kinetics.product_stoich_coeffs(gas) # Stoichiometric coefficient matrix {Products}
v_r = ct.Kinetics.reactant_stoich_coeffs(gas) # Stoichiometric coefficient matrix {Reactants}
vk = v_p - v_r # Net stoichiometric coefficients [v_p - v_r]
Ru = ct.gas_constant/1000 # Universal gas constant [J/mol.K]

# Set the ODE solver
IC = np.array([X0,T0])
sol = odeint(RMGCanODE,IC,t0,args=(gas,p0,patm,NM),atol=1e-12,rtol=1e-3)


gas()
sol

У меня есть аналогичный код в MATLAB, который отлично работает. Итак, я считаю, что метод правильный, то есть химическая и термодинамическая части. Я пробовал несколько способов упаковать начальные условия и возвращаемую часть функции ODE. Однако самая последняя итерация приводит к следующей ошибке TypeError: невозможно преобразовать данные массива из dtype('O') в dtype('float64') в соответствии с правилом "safe".

Где я мог ошибиться?

-ps- Чтобы решить эту проблему в вашей системе, вам необходимо установить Cantera (инструкции по установке на https://cantera.org/install/index.html).