# 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. – Lutz Lehmann Oct 30 '20 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. – Wolfgang Bangerth Oct 31 '20 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? – user606273 Nov 2 '20 at 14:51

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.