решающая система связанных од с одеинт

Question

Я использую систему ода для моделирования обжарки кофейных зерен для задания класса.Уравнения приведены ниже.

Все параметры (кроме 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()

У вас есть идеи, почему шаг интеграции не работает?

Спасибо !!

Answer 1

Вы неоднократно решаете систему уравнений только для одного момента времени.

Из документации odeint команда odeint принимает аргумент t, который:

Последовательность моментов времени, для которых нужно решить для y. Начальное значение точки должно быть первым элементом этой последовательности.

Поскольку вы передаете [t+dt] в odeint, существует только один момент времени, поэтому вы получаете только одно значение, которое просто является вашим начальным условием.

Правильный способ использования odeint похож на следующее:

output = odeint(deriv, [T_b, X_b], np.linspace(0,600,600))

Здесь output, опять же согласно документации это:

Массив, содержащий значение y для каждого требуемого времени в t, с начальным значением y0 в первой строке.