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$ – Lutz Lehmann Oct 30 '20 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$ – Wolfgang Bangerth Oct 31 '20 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 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.

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, privacy policy and cookie policy

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