0

I have been trying to solve the following system of three coupled differential equations using odeint and solve_ivp. (This is not the exact system. I have modified it just to improve the readability of the question.)

def model(x,t):    # x is a list containing the values of u,v,w as they change with time t 
    u = x[0]
    v = x[1]
    w = x[2]
    
    dudt = u - u*v                           
    dvdt = -v + u*v - v*w**2          
    dwdt = -w + 2*v*w
    
    return [dudt, dvdt, dwdt]

This works fine. But now I want to modify this in the following way: whenever either of u,v or w goes below a threshold, it is reset to zero, and then let the system evolve. This should happen every time any of these three goes before the threshold automatically. The 'rules' of evolving the system remain the same.

I have tried modifying the code by doing this:

def model(x,t):    # x is a list containing the values of u,v,w as they change with time t 
    u = x[0]
    v = x[1]
    w = x[2]
     
    if u < u_threshold:
        x[0] = 0
    
    dudt = u - u*v                           
    dvdt = -v + u*v - vw**2          
    dwdt = -w + 2*v*w
    
    return [dudt, dvdt, dwdt]

I have shown it only for u, but you get the idea. This does not seem to work.

Please note that I cannot afford to stop the simulation every time any variable hits the threshold value, as this is only a toy-model. Later on, I will generalise this to systems of hundreds of coupled equations.

Singh
  • 1
  • 1
  • What do you mean with "stop the simulation"? The simulation is a sequence of time steps, one could say it "temporarily stops" in every step. – Lutz Lehmann Sep 28 '20 at 07:23
  • What is the reasoning behind these jumps? Is that to avoid some numerical complications or is this a hard consequence of the theory behind the simulation? – Lutz Lehmann Sep 28 '20 at 07:31
  • By "stop the simulation", I meant permanently stopping it. – Singh Sep 28 '20 at 08:08
  • The reasoning: This system is modelling a biological system, where u, v and w are the density of different species of organisms in the ecosystem. There is an interesting effect seen in biological systems: if the density goes below a threshold value, usually the species goes extinct. – Singh Sep 28 '20 at 08:10
  • Ok, that makes some sense, you want to prevent that a species recovers from a population of 0.1 individuals (dividing the in-dividuum, the non-divisible). The regular way would be to run the simulation for a fixed time, a month or a year, and then prune the state vector and restart with the modified state. This is not noticeably slower than a continuous simulation, array operations exist to concatenate the resulting arrays of states. See https://stackoverflow.com/a/59054871/3088138 for an implementation of this for a slightly different reason. – Lutz Lehmann Sep 28 '20 at 08:18
  • It does not matter here, but just to clarify- species go extinct even when there are multiple individuals. For example, if you have too few lions distributed in a large forest, the probability of two individuals mating can be too low. – Singh Sep 28 '20 at 08:52

1 Answers1

0

It does not work, as the state-vector x is read-only, there is no side-effect transmitting changes back to the state vector of the integrator.

And even if it did work, it would be a bad idea:

  • The ODE function model is mostly not called on points of the solution curve, and then in the higher order RK methods also not with increasing times.

  • Such a jump would make the ODE highly non-smooth destroying all the assumptions that are used in the time-step-size controller. The result of this can be varied, but all non-favorable towards a minimum-effort integration.

So yes, the best course is indeed to stop the simulation using the event mechanism of solve_ivp or to fashion your own time loop based on the stepper classes that are also behind solve_ivp.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51