python中的函数和odeint的问题

Question

几个月来，我开始使用python，考虑到它的巨大优势。但最近，我从scipy使用odeint来求解一个微分方程组。但在整合过程中，实现的功能无法按预期工作。
在这种情况下，我想要求解一个微分方程组，其中一个初始条件（x [0]）在4-5之间变化，取决于变量在积分过程中达到的值（它被编程通过if结构的函数内部）。

#Control of oxygen 
    SO2_lower=4 
    SO2_upper=5 
    if x[0]<=SO2_lower: 
     x[0]=SO2_upper 

当函数由odeint，在该函数内的一些代码行可以被排除，即使当函数改变x的值[0]。这里是我所有的代码：

import numpy as np 
    from scipy.integrate import odeint 
    import matplotlib.pyplot as plt 
    plt.ion() 
    # Stoichiometric parameters 
    YSB_OHO_Ox=0.67       #Yield for XOHO growth per SB (Aerobic) 
    YSB_Stor_Ox=0.85       #Yield for XOHO,Stor formation per SB (Aerobic) 
    YStor_OHO_Ox=0.63       #Yield for XOHO growth per XOHO,Stor (Aerobic) 
    fXU_Bio_lys=0.2       #Fraction of XU generated in biomass decay 
    iN_XU=0.02        #N content of XU 
    iN_XBio=0.07        #N content of XBio 
    iN_SB=0.03        #N content of SB 
    fSTO=0.67         #Stored fraction of SB 

    #Kinetic parameters 
    qSB_Stor=5        #Rate constant for XOHO,Stor storage of SB 
    uOHO_Max=2        #Maximum growth rate of XOHO 
    KSB_OHO=2         #Half-saturation coefficient for SB 
    KStor_OHO=1        #Half-saturation coefficient for XOHO,Stor/XOHO 
    mOHO_Ox=0.2        #Endogenous respiration rate of XOHO (Aerobic) 
    mStor_Ox=0.2        #Endogenous respiration rate of XOHO,Stor (Aerobic) 
    KO2_OHO=0.2        #Half-saturation coefficient for SO2 
    KNHx_OHO=0.01        #Half-saturation coefficient for SNHx 

    #Other parameters 
    DT=1/86400.0 

    def f(x,t): 
     #Control of oxygen 
     SO2_lower=4 
     SO2_upper=5 
     if x[0]<=SO2_lower: 
      x[0]=SO2_upper 
     M=np.matrix([[-(1.0-YSB_Stor_Ox),-1,iN_SB,0,0,YSB_Stor_Ox], 
      [-(1.0-YSB_OHO_Ox)/YSB_OHO_Ox,-1/YSB_OHO_Ox,iN_SB/YSB_OHO_Ox-iN_XBio,0,1,0], 
      [-(1.0-YStor_OHO_Ox)/YStor_OHO_Ox,0,-iN_XBio,0,1,-1/YStor_OHO_Ox], 
      [-(1.0-fXU_Bio_lys),0,iN_XBio-fXU_Bio_lys*iN_XU,fXU_Bio_lys,-1,0], 
      [-1,0,0,0,0,-1]]) 
     R=np.matrix([[DT*fSTO*qSB_Stor*(x[0]/(KO2_OHO+x[0]))*(x[1]/(KSB_OHO+x[1]))*x[4]], 
      [DT*(1-fSTO)*uOHO_Max*(x[0]/(KO2_OHO+x[0]))*(x[1]/(KSB_OHO+x[1]))* (x[2]/(KNHx_OHO+x[2]))*x[4]], 
      [DT*uOHO_Max*(x[0]/(KO2_OHO+x[0]))*(x[2]/(KNHx_OHO+x[2]))*((x[5]/x[4])/(KStor_OHO+(x[5]/x[4])))*(KSB_OHO/(KSB_OHO+x[1]))*x[4]], 
      [DT*mOHO_Ox*(x[0]/(KO2_OHO+x[0]))*x[4]], 
      [DT*mStor_Ox*(x[0]/(KO2_OHO+x[0]))*x[5]]]) 

     Mt=M.transpose() 
     MxRm=Mt*R 
     MxR=MxRm.tolist() 
     return ([MxR[0][0], 
       MxR[1][0], 
       MxR[2][0], 
       MxR[3][0], 
       MxR[4][0], 
       MxR[5][0]]) 
    #ODE solution 
    t=np.linspace(0.0,3600,3600) 
    #Initial conditions 
    y0=np.array([5,176,5,30,100,5]) 
    Var=odeint(f,y0,t,args=(),h0=1,hmin=1,hmax=1,atol=1e-5,rtol=1e-5) 
    Sol=Var.tolist() 
    plt.plot(t,Var[:,0]) 

非常感谢您提前致谢!!!!!

Answer 1

简短回答：

您不应该修改ODE函数内的输入状态向量。相反，请尝试以下和验证结果：

x0 = x[0] 
if x0<=SO2_lower: 
    x0=SO2_upper 
# use x0 instead of x[0] in the rest of this function body 

我想，这是你的问题，但我不知道，因为你没有解释究竟是错的结果。而且，你不会改变“初始条件”。初始条件是

y0=np.array([5,176,5,30,100,5]) 

你只是改变输入状态向量。

详细的解答：

你odeint积分器可能使用的高阶自适应龙格 - 库塔方法之一。该算法需要多个ODE函数评估来计算单个积分步骤，因此改变输入状态向量可能导致未定义的结果。在C++ boost :: odeint中，这甚至不可能这样做，因为输入变量是“const”。然而，Python并不像C++那样严格，我认为有可能无意中造成这种错误（尽管我没有尝试过）。

编辑：

好的，我明白你想达到什么。

你的变量x [0]是通过模块化代数约束，并且因此不可能在形式

x' = f(x,t) 

这是常微分方程的可能的定义中的一个来表示，即ondeint库是指解决。然而，在这里可以使用很少的“黑客”来绕过这个限制。

一种可能性是使用固定的步骤和低位（因为你需要知道的高阶解算器，你实际上是在算法的一部分，看到RK4例如）求解器和改变你的DX [0]方程（在你的代码是MXR [0] [0]元素）到：

# at the beginning of your system 
if (x[0] > S02_lower): # everything is normal here 
    x0 = x[0] 
    dx0 = # normal equation for dx0 
else: # x[0] is too low, we must somehow force it to become S02_upper again 
    dx0 = (x[0] - S02_upper)/step_size // assuming that x0_{n+1} = x0_{n} + dx0*step_size 
    x0 = S02_upper 
# remember to use x0 in the rest of your code and also remember to return dx0 

不过，我不推荐这种技术，因为它使你强烈地依赖于算法，你必须知道确切的步长（尽管，我可能会建议它为饱和约束）。另一种可能性是在同一时间执行单积分步，每次纠正你X0是必要的：

// 1 do_step(sys, in, dxdtin, t, out, dt); 
// 2 do something with output 
// 3 in = out 
// 4 return to 1 or finish 

对不起，C++的语法，这里是详尽的文档（​​），这里是等值的蟒蛇（Python ode class）。 C++ odeint接口更适合你的任务，但是你可能在python中完全一样。只需寻找：

integrate(t[, step, relax]) 
set_initial_value(y[, t]) 

in docs。