Problems With A Function And Odeint In Python

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 dur

Solution 1:

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 for example) solver and change your dx[0] equation (in your code it is MxR[0][0] element) to:

# at the beginning of your systemif (x[0] > S02_lower): # everything is normal here
    x0 = x[0]
    dx0 = # normal equation for dx0else: # 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), and here is its equivalent in python (Python ode class). 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.