Problems with a function and odeint in python
Question: For a few months I started working with python, considering the great
advantages it has. But recently, i used odeint from scipy to solve a system of
differential equations. But during the integration process the implemented
function doesn't work as expected.  
In this case, I want to solve a system of differential equations where one of
the initial conditions (x[0]) varies (between 4-5) depending on the value that
the variable reaches during the integration process (It is programmed inside
of the function by means of the if structure).

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

When the function is used by odeint, some lines of code inside the function
are obviated, even when the functions changes the value of x[0]. Here is all
my code:

    
    
        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]) 
    

Thanks very much in advance!!!!!


Answer: Short answer:

You should not modify input state vector inside your ODE function. Instead try
the following and verify your results:

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

I suppose that this is your problem, but I am not sure, since you did not
explain what exactly was wrong with the results. Moreover, you do not change
"initial condition". Initial condition is

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

you just change the input state vector.

Detailed answer:

Your odeint integrator is probably using one of the higher order adaptive
Runge-Kutta methods. This algorithm requires multiple ODE function evaluations
to calculate single integration step, therefore changing the input state
vector may lead to undefined results. In C++ boost::odeint this is even not
possible to do so, because input variable is "const". Python however is not as
strict as C++ and I suppose that it is possible to make this kind of bug
unintentionally (I did not try it, though).

EDIT:

OK, I understand what you want to achieve.

Your variable x[0] is constrained by modular algebra and it is not possible to
express in the form

    
    
    x' = f(x,t)
    

which is one of the possible definitions of the Ordinary Differential
Equation, that ondeint library is meant to solve. However, few possible
"hacks" can be used here to bypass this limitation.

One possibility is to use a fixed step and low order (because for higher order
solvers you need to know, which part of the algorithm you are actually in, see
[RK4](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge.E2.80.93Kutta_method)
for example) solver and change your dx[0] equation (in your code it is
MxR[0][0] element) to:

    
    
    # 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
    

However, I do not recommend this technique, because it makes you strongly
dependent on the algorithm and you must know the exact step size (although, I
may recommend it for saturation constraints). Another possibility is to
perform a single integration step at a time and correct your x0 each time it
is necessary:

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

Sorry for C++ syntax, here is the exhaustive documentation ([C++ odeint
steppers](http://www.boost.org/doc/libs/master/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html)),
and here is its equivalent in python ([Python ode
class](http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode)).
C++ odeint interface is better for your task, however you may achieve exactly
the same in python. Just look for:

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

in docs.

