ScipyでODEの集合を解く

Question

システムに到達するまでの時間の経過とともに変化する細胞、基質および生成物（すなわち、）の濃度を予測するアルゴリズム（scipy.integrate.odeint（）を使用）を開発しようとしています。定常状態（約100時間または200時間）。バイオリアクター中の細胞の初期濃度は0.1 /であり、最初は反応器中にグルコースまたは生成物が存在しない。 0.01 /ℎと0.25/betweenの間のさまざまな流量のアルゴリズムをテストし、流量が製品生産に及ぼす影響を分析したいと考えています（⋅in /ℎ）。最終的には、製品生産率（y軸）対流量（y軸）をx軸に示すプロットを生成したいと考えています。私の目標は、最大（または臨界）生産率をもたらす流速を推定することです。私は値を取得するには、この修正を試み、現在の状態がPの定常状態値を生成するための方法であることがあるようには思えないプログラムで

from scipy.integrate import odeint 
import numpy as np 

# Constants 
u_max = 0.65 
K_s = 0.14 
K_1 = 0.48 
V = 2 
X_in = 0 
S_in = 4 
Y_s = 0.38 
Y_p = 0.2 

# Variables 
# Q - Flow Rate (L/h), value between 0.01 and 0.25 that produces best Q * P 
# X - Cell Concentration (g/L) 
# S - The glucose concentration (g/L) 
# P - Product Concentration (g/L) 

# Equations 
def func_dX_dt(X, t, S): 
    u = (u_max)/(1 + (K_s/S)) 
    dX_dt = (((Q * S_in) - (Q * S))/V) + (u * X) 
    return dX_dt 

def func_dS_dt(S, t, X): 
    u = (u_max)/(1 + (K_s/S)) 
    dS_dt = (((Q * S_in) - (Q * S))/V) - (u * (X/Y_s)) 
    return dS_dt 

def func_dP_dt(P, t, X, S): 
    u = (u_max)/(1 + (K_s/S)) 
    dP_dt = ((-Q * P)/V) - (u * (X/Y_p)) 
    return dP_dt 

t = np.linspace(0, 200, 200) 

# Q placeholder 
Q = 0.01 

# Attempt to solve the Ordinary differential equations 
sol_dX_dt = odeint(func_dX_dt, 0.1, t, args=(S,)) 
sol_dS_dt = odeint(func_dS_dt, 0.1, t, args=(X,)) 
sol_dP_dt = odeint(func_dP_dt, 0.1, t, args=(X,S)) 

：これは、これまでの私のコードですX.

sol_dX_dt = odeint(func_dX_dt, 0.1, t, args=(odeint(func_dS_dt, 0.1, t, args=(X,)),)) 

のそれはエラーを生成します。

THIで

NameError: name 'X' is not defined 

私は前進する方法がわかりません。

（編集1：追加オリジナル方程式）

まず式

第2式と第三式

Answer 1

あなたはに関数を適用する必要はありません各部分は、以下に示すように、派生物のタプルを返します：

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

Q = 0.01 
V = 2 
Ys = 0.38 
Sin = 4 
Yp = 0.2 
Xin = 0 
umax = 0.65 
Ks = 0.14 
K1 = 0.48 

def mu(S, umax, Ks, K1): 
    return umax/((1+Ks/S)*(1+S/K1)) 

def dxdt(x, t, *args): 
    X, S, P = x 
    Q, V, Xin, Ys, Sin, Yp, umax, Ks, K1 = args 
    m = mu(S, umax, Ks, K1) 
    dXdt = (Q*Xin - Q*X)/V + m*X 
    dSdt = (Q*Sin - Q*S)/V - m*X/Ys 
    dPdt = -Q*P/V - m*X/Yp 
    return dXdt, dSdt, dPdt 

t = np.linspace(0, 200, 200) 
X0 = 0.1 
S0 = 0.1 
P0 = 0.1 
x0 = X0, S0, P0 
sol = odeint(dxdt, x0, t, args=(Q, V, Xin, Ys, Sin, Yp, umax, Ks, K1)) 

plt.plot(t, sol[:, 0], 'r', label='X(t)') 
plt.plot(t, sol[:, 1], 'g', label='S(t)') 
plt.plot(t, sol[:, 2], 'b', label='P(t)') 
plt.legend(loc='best') 
plt.xlabel('t') 
plt.grid() 
plt.show() 

出力：あなたはこの内のコードので、方程式（多分式の画像がこの場合は正しい）の初期条件、および他の式を置くことができ