Проблемы с функцией и odeint в Python

В течение нескольких месяцев я начал работать с python, учитывая его огромные преимущества. Но недавно я использовал odeint от scipy для решения системы дифференциальных уравнений. Но в процессе интегрирования реализованная функция работает не так, как ожидалось.
В этом случае я хочу решить систему дифференциальных уравнений, в которой одно из начальных условий (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])

Большое спасибо заранее!!!!!


python-3.x odeint
person Yusmel González Hernández    schedule 22.01.2016    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Короткий ответ:

Вы не должны изменять вектор состояния ввода внутри функции 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])

вы просто меняете вектор состояния ввода.

Подробный ответ:

Ваш единственный интегратор, вероятно, использует один из адаптивных методов Рунге-Кутты более высокого порядка. Этот алгоритм требует множественных вычислений функции 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). Интерфейс C ++ odeint лучше подходит для вашей задачи, однако вы можете добиться того же в python. Просто ищите:

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

в документах.

person pptaszni    schedule 02.02.2016
comment
Большое спасибо за ваш ответ. В моем случае переменная x [0] изменяется от 4 до 5, и поэтому каждый раз, когда x [0] достигает значения, равного или меньшего 4, его необходимо снова сбрасывать. Есть какие-то варианты решения этой проблемы? - person Yusmel González Hernández; 04.02.2016
comment
Это вполне возможно, но требует некоторых небольших хитростей. Пожалуйста, прочтите отредактированный пост выше. - person pptaszni; 04.02.2016
comment
Я постараюсь использовать ODE, как вы предлагаете. Большое спасибо, мой друг. - person Yusmel González Hernández; 05.02.2016