如何使用 scipy.odeint 求解微分方程组

Question

我想用 odeint 求解方程组，但出现以下错误：

  File "C:", line 45, in <module>
    C_B = odeint(dC_Bdt,C_B0,t)

  File "C:\Anaconda3\envs\ChemEng\lib\site-packages\scipy\integrate\odepack.py", line 233, in odeint
    int(bool(tfirst)))

RuntimeError: The array return by func must be one-dimensional, but got ndim=2.

我的代码是：

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

# Create time domain
t = np.linspace(0,100,100)

# Parameters 
A = 1*10**(13) # Arrhenius constant
T = 293.15 # Temperature [K]
E_a= 80000 # Activation energy [J/mol]
R = 8.31 # Ideal gas constant [J/molK]
rho = 1000 # density [kg/m3]
F_in = 0.2 # Inlet flowrate [m3/s]
h = 2.1 # Height of reactor
A_= 1 # Cross-sectional area of reactor [m2]


# Initial condition
C_A0 = 3 # Initial concentration [mol/m3]
C_B0 = 0 # Initial concentration [mol/m3]
m0 = 0 # Initial mass in tank [kg]

def dmdt(F_out,t): # Mass balance
    return rho*(F_in-F_out)

def dC_Adt(C_A,t): # Concentration balance for A
    dC_Adt = (F_in*C_A-F_out*C_A)/V-k*C_A
    return dC_Adt

def dC_Bdt(C_B,t): # Concentration balance for B
    dC_Bdt = (F_in*C_B-F_out*C_B)/V+k*C_A
    return dC_Bdt #<-- needs to be 1D but is 2D


V = A_*h # Reactor volume [m3]
k=A*np.exp(-E_a/(R*T)) # Reaction rate constant
F_out = F_in # Steady state
dmdt = odeint(dmdt,m0,t)
C_A = odeint(dC_Adt,C_A0,t)
C_B = odeint(dC_Bdt,C_B0,t)

# Plot
plt.figure()
plt.plot(t,C_A,'b-',label='C_A')
plt.plot(t,C_B,'r--',label='C_B')
plt.legend(loc='best')
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Concentration [mol/m3]')

在以下位置提出了类似的问题：How to fix the error: The array return by func must be one-dimensional, but got ndim=2

但答案是更改函数依赖项的顺序。那对我不起作用。

我认为问题出在我使用不同函数 (C_A) 的输出来计算 C_B。那么 dC_Bdt 是否会以某种方式将 C_A 和 C_B 作为输出，或者至少以某种方式将它们链接起来，这就是为什么在使用 odeint 时出现该错误？

不确定从这里开始做什么？ :( 非常感谢

Answer 1

我修改了您的代码，以便求解三个耦合方程：

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

# Parameters 
A = 1e13 # Arrhenius constant
T = 293.15 # Temperature [K]
E_a = 80000 # Activation energy [J/mol]
R = 8.31 # Ideal gas constant [J/molK]
rho = 1000 # density [kg/m3]
F_in = 0.2 # Inlet flowrate [m3/s]
h = 2.1 # Height of reactor
A_ = 1 # Cross-sectional area of reactor [m2]

V = A_*h  # Reactor volume [m3]
k = A*np.exp(-E_a/(R*T)) # Reaction rate constant
F_out = F_in # Steady state

def dUdt(U, t):
    m, C_A, C_B = U

    dmdt = rho*(F_in - F_out)
    dC_Adt = ( F_in*C_A - F_out*C_A )/ V-k*C_A
    dC_Bdt = ( F_in*C_B - F_out*C_B )/ V+k*C_A

    return [dmdt, dC_Adt, dC_Bdt]

# Create time domain
t_span = np.linspace(0, 100, 30)

# Initial condition
C_A0 = 3 # Initial concentration [mol/m3]
C_B0 = 0 # Initial concentration [mol/m3]
m0 = 0 # Initial mass in tank [kg]

Uzero = [m0, C_A0, C_B0]

solution = odeint(dUdt, Uzero, t_span)

# plot
plt.plot(t_span, solution[:, 0], label='masse');
plt.plot(t_span, solution[:, 1], label='C_A');
plt.plot(t_span, solution[:, 2], label='C_B');
plt.legend();
plt.xlabel('time');

如何使用 scipy.odeint 求解微分方程组

How to solve a system of differential equations using scipy.odeint

python

scipy

odeint