CODEWITHSUNDEEP
pythonpython-3.xscipydifferential-equationsodeint
Oniow
Oniow

Reputation: 329

Controlling output limits in Scipy's ODEINT

I have the following simplified pseudocode setup which involves a system of ODEs, which I attempt to solve with scipy's odeint.

from scipy.integrate import odeint

def diff_func(y, time, parms):

    # Do stuff with parms that depends upon y and t.
    new_parms = other_funcs(y,time, parms)

    # Now calculate derivatives
    dy1_dt = dy1_func(y, new_parms)
    dy2_dt = dy2_func(y, new_parms)

# Setup up initial conditions
y_0 = [y_1_0, y_2_0]
time = np.linspace(0, 1000, 1000)
parms = list_o_constants

# Solve diffeq
yout, info = odeint(diff_func, y_0, time, args=(parms,), full_output=True)

Depending upon the inputs I give it, I can lead to nonsensical results - or crashes. I have narrowed down the problem to the y inputs go negative occasionally, since one of them is temperature - it leads to all kinds of issues. The reason why it is going negative in the first place is unclear, and might be a stability issue.

I have attempted to change the diff function to include a check like:

def diff_func(y, time, parms):
    if y[1] <= 0.1:
        y[1] = 0.1

This seems to work okay, I am still testing it. But, I wanted to ask if: this is going to cause problems down the line, and if there is a better way to force odeint to only allow positive outputs (i.e., change step size if a negative is found).

Thanks!

Upvotes: 3

Views: 2102

Answers (1)

pptaszni
pptaszni

Reputation: 8217

Please first look at my other answer that I posted recently.

What you want do do is to add a saturation constraint to your ODE. There are a few possibilities to achieve this behavior, however I recommend this one:

# calculate dy[1]
y1 = y[1]
if dy[1] < 0 and y1 <= 0.1 # or 0.0, depends on accuracy you wish to maintain
    y1 = 0
    dy1 = 0 # this prevents your variable from going further down
else
    # normal routine

Let me know, if it suits you.

Best regards.

Upvotes: 1

Related Questions