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

本文介绍了如何使用 scipy.odeint 求解微分方程组的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述

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

I want to solve a system of equations using odeint and I get the following error:

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]')

一个类似的问题被问到:如何修复错误:func返回的数组必须是一维的，但得到了ndim=2

A similar question was asked at: How to fix the error: The array return by func must be one-dimensional, but got ndim=2

但答案是改变函数依赖的顺序.那对我不起作用.

But the answer to that was to change the order of the function's dependencies. That doesn't work for me.

我认为问题源于我使用不同函数 (C_A) 的输出来评估 C_B.那么可能是 dC_Bdt 以某种方式将 C_A 和 C_B 作为输出，或者至少以某种方式将它们链接起来，这就是为什么在使用 odeint 时出现错误?

I think the problem arises from the fact that I use an output of a different function (C_A), to evaluate C_B. So could it be that dC_Bdt somehow has C_A and C_B as ouputs or at least links them somehow which is why when using odeint I get that error?

不知道从这里做什么?:(非常感谢

Not sure what to do from here? :( Many Thanks

推荐答案

我修改了你的代码，为了解决三个耦合方程:

I modified your code, in order to solve the three coupled equations:

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');