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):
 x, u= z
 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))
  • $\begingroup$ 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. $\endgroup$ Commented Oct 30, 2020 at 13:59
  • $\begingroup$ 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. $\endgroup$ Commented Oct 31, 2020 at 0:36
  • $\begingroup$ 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? $\endgroup$
    – user606273
    Commented Nov 2, 2020 at 14:51
  • $\begingroup$ 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? $\endgroup$
    – R. Engel
    Commented May 23, 2022 at 8:26

1 Answer 1


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.

Finally, you can hack your way around by simply projecting after each step the solution components back on the boundaries if they are outside of the allowed interval. This is a pragmatic solution, which may not always be well-behaved...

One important thing is that your solution limitation must not "fight" the ODE, i.e. you should not force one variable u to be >= 0 if f(u=0) is < 0 for example, as this will only cause trouble (oscillations, unnecessarily low time step values) in the integration. One solution is to modify your ODE function so that the time derivatives of the variables slowly decay to 0 as they are getting closer to their bounds.

You can read this paper to learn more.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.