用 odeint 求解耦合 odes 系统
solving system of coupled odes with odeint
我正在使用 ode 系统为 class 作业模拟咖啡豆烘焙。等式如下。
参数(X_b和T_b除外)均为常量。
当我尝试使用 odeint 来解决这个系统时,它给出了一个常量 T_b 和 X_b 配置文件(这在概念上没有意义)。
下面是我使用的代码
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
# Write function for bean temperature T_b differential equation
def deriv(X,t):
T_b, X_b = X
dX_b = (-4.32*10**9*X_b**2)/(l_b**2)*np.exp(-9889/T_b)
dT_b = ((h_gb*A_gb*(T_gi - T_b))+(m_b*A_arh*np.exp(-H_a/R_g/T_b))+
(m_b*lam*dX_b))/(m_b*(1.099+0.0070*T_b+5*X_b)*1000)
return [dT_b, dX_b]
# Establish initial conditions
t = 0 #seconds
T_b = 298 # degrees K
X_b = 0.1 # mass fraction of moisture
# Set time step
dt = 1 # second
# Establish location to store data
history = [[t,T_b, X_b]]
# Use odeint to solve DE
while t < 600:
T_b, X_b = odeint(deriv, [T_b, X_b], [t+dt])[-1]
t += dt
history.append([t,T_b, X_b])
# Plot Results
def plot_history(history, labels):
"""Plots a simulation history."""
history = np.array(history)
t = history[:,0]
n = len(labels) - 1
plt.figure(figsize=(8,1.95*n))
for k in range(0,n):
plt.subplot(n, 1, k+1)
plt.plot(t, history[:,k+1])
plt.title(labels[k+1])
plt.xlabel(labels[0])
plt.grid()
plt.tight_layout()
plot_history(history, ['t (s)','Bean Temperature $T_b$ (K)', 'Bean Moisture Content $X_b$'])
plt.show()
您知道为什么集成步骤不起作用吗?
谢谢!!
您只是针对一个时间点重复求解方程组。
从 odeint documentation 开始,odeint
命令接受一个参数 t
,即:
A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.
由于您将 [t+dt]
传递给 odeint
,因此只有一个时间点,因此您只能返回一个值,这只是您的初始条件。
正确的使用方法odeint
类似如下:
output = odeint(deriv, [T_b, X_b], np.linspace(0,600,600))
这里output
,再次根据文档是:
Array containing the value of y for each desired time in t, with the initial value y0 in the first row.
我正在使用 ode 系统为 class 作业模拟咖啡豆烘焙。等式如下。
参数(X_b和T_b除外)均为常量。
当我尝试使用 odeint 来解决这个系统时,它给出了一个常量 T_b 和 X_b 配置文件(这在概念上没有意义)。
下面是我使用的代码
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
# Write function for bean temperature T_b differential equation
def deriv(X,t):
T_b, X_b = X
dX_b = (-4.32*10**9*X_b**2)/(l_b**2)*np.exp(-9889/T_b)
dT_b = ((h_gb*A_gb*(T_gi - T_b))+(m_b*A_arh*np.exp(-H_a/R_g/T_b))+
(m_b*lam*dX_b))/(m_b*(1.099+0.0070*T_b+5*X_b)*1000)
return [dT_b, dX_b]
# Establish initial conditions
t = 0 #seconds
T_b = 298 # degrees K
X_b = 0.1 # mass fraction of moisture
# Set time step
dt = 1 # second
# Establish location to store data
history = [[t,T_b, X_b]]
# Use odeint to solve DE
while t < 600:
T_b, X_b = odeint(deriv, [T_b, X_b], [t+dt])[-1]
t += dt
history.append([t,T_b, X_b])
# Plot Results
def plot_history(history, labels):
"""Plots a simulation history."""
history = np.array(history)
t = history[:,0]
n = len(labels) - 1
plt.figure(figsize=(8,1.95*n))
for k in range(0,n):
plt.subplot(n, 1, k+1)
plt.plot(t, history[:,k+1])
plt.title(labels[k+1])
plt.xlabel(labels[0])
plt.grid()
plt.tight_layout()
plot_history(history, ['t (s)','Bean Temperature $T_b$ (K)', 'Bean Moisture Content $X_b$'])
plt.show()
您知道为什么集成步骤不起作用吗?
谢谢!!
您只是针对一个时间点重复求解方程组。
从 odeint documentation 开始,odeint
命令接受一个参数 t
,即:
A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.
由于您将 [t+dt]
传递给 odeint
,因此只有一个时间点,因此您只能返回一个值,这只是您的初始条件。
正确的使用方法odeint
类似如下:
output = odeint(deriv, [T_b, X_b], np.linspace(0,600,600))
这里output
,再次根据文档是:
Array containing the value of y for each desired time in t, with the initial value y0 in the first row.