python中函数和odeint的问题-Java 学习之路

几个月后，我开始使用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])

首先十分感谢！！！！！

1 回答

0
简短回答：

您不应该修改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集成商可能正在使用其中一种更高阶的自适应Runge-Kutta方法 . 该算法需要多个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语法，这里是详尽的文档（C++ odeint steppers），这是它在python（Python ode class）中的等价物 . C odeint接口更适合您的任务，但是您可以在python中实现完全相同 . 只需寻找：
```
integrate(t[, step, relax])
set_initial_value(y[, t])
```
在文档中 .
回复于 2024-04-30T14:55:38+08:00

python中函数和odeint的问题

1 回答

相关问题