# Is there a way to bound the values of a variable when using scipy.integrate.solve_ivp in python?

I want to solve an IVP in python with two variables, x and u, but I need the values of u to be between 0 and 1. Right now it is giving me a solution with negative values for u. Here is the code I have.

def model_case_3(t, z, Kmax, k, b, list_Kmax, sigma):
a=1.5
x, u= z
m=1
r = list_Kmax[0][0]
dxdt =  (x)*(r*(1-a*u**2)*(1-x/(Kmax*(1-0.999*u**2)))-m/(k+b*u)-0.05)
dudt = sigma*((-2*a*(b**2)*r*(u**3)+4*a*b*k*r*(u**2)+2*a*(k**2)*r*u-b*m)/((b*u+k)**2))
return [dxdt, dudt]

sol = solve_ivp(fun=model_case_3, t_span=[scaled_days[i][j][0], scaled_days[i][j][-1]], y0=[scaled_pop[i][j][0], list_u[0][0][0]], t_eval=scaled_days[i][j],  args=(list_Kmax[0][0], k0, b0, list_b, sigma0))

• You would need to change the ODE so that the vector field does not lead into the region of negative $u$. That can happen internally by modifying the function in some way that is not too singular, or by defining an event and then switching to a different ODE function. Commented Oct 30, 2020 at 13:59
• What @LutzLehmann is saying, is that the ODE is what it is. It's going to produce a solution that is what it is. If you want something else, you need a different ODE. Commented Oct 31, 2020 at 0:36
• Yes, I see your point, but is there any way that I can kind of "proyect" the values of u so that if they are bigger than 1 they are taken to be 1 or something like this? Commented Nov 2, 2020 at 14:51
• I have the same problem, but my system has no meaningful solution outside certain parameter bounds so that when the solver goes there I often have to throw an error. Yes, small tolerances can solve this but become computationally prohibitive. Is there no way of passing constrains to the solver? Commented May 23, 2022 at 8:26

If the exact solution indeed stays within [0,1], the solver may still resolve too coarsely the dynamics and "jump" over the physical bounds. One way to solve this is to use lower absolute and relative tolerances in your solve_ivp call: solve_ivp(..., atol=1e-9, rtol=1e-9, ... for example.