1

I need to solve a system of differential equations and I have a function that returns dx/dt and du/dt, but I need to add the constraint that u is in the interval [0,1]. How can I do that? So far I've been using odeint from scipy.integrate to solve it. This is the model:

def model_case_1(z, t, Kmax):
    '''The input z corresponds to the current state of the system, z = [x, u, m] 
    (x is the population, u is the resistance rate and m whether there is treatment or not). 
    t is the current time.
    Kmax corresponds to the unknown parameter.
    '''
    
    x, u= z
    if t>0:
      m=1
    else:
      m=0
    dxdt =  x*(r*(1-x//(Kmax*(1-u)))-m//(k+b*u)-d)
    dudt = sigma*(m*b//((k+b*u)**2)-r*x//(Kmax*(1-u)**2))
    return [dxdt, dudt]

Then I estimate Kmax and I solve the ODE with:

    y = odeint(model_case_1, [x0, u0], teval, args=(Kmax0,))

But then, of course, I get values for u bigger than 1, which cannot happen.

user606273
  • 143
  • 6
  • Is the use of integer division intentional? Is the restriction of `u` in any way natural in the theoretical model? – Lutz Lehmann Sep 10 '20 at 21:33
  • Thanks for pointing that out! I just changed it because I was getting an error and then forgot about it, but yeah, it should be normal division – user606273 Sep 11 '20 at 10:23

0 Answers0