I have been trying to solve a system of differential equations in Julia. I have run into a few problems and am not clear on what I am doing wrong.
The system I am working on is a model of the immune system. This is represented by 9 differential equations. One of the parameters is time-dependent and is a step function representing the lifetime of an infection in the immune system. This makes the system discontinuous at the time of entry and removal of the infection. To address this, I used the options tstops and callbacks in the solver. This system seems to be stiff and so I used the Rosenbrock23 algorithm to solve this.
Because I do not know the parameters, they are assigned random values. For independent simulations (for different parameter sets), I sometimes run into negative solutions. I tried using the callback PositiveDomain but that does not help. Even for simulations where Rosenbrock23 algorithm handles the integration quite well, adding the callback PositiveDomain gives negative solution and aborts with the warning:
WARNING: dt <= dtmin. Aborting. If you would like to force continuation with dt=dtmin, set force_dtmin=true
I understand that the step size can’t be reduced anymore and that’s why it’s aborting but I don’t understand why adding the callback PositiveDomain would do this. Or is my implementation faulty?
In cases, where I do get negative values for some variables, isoutofbound solver option also doesn’t seem to make any difference.
What am I doing wrong here?
The code: https://www.pastiebin.com/5b0560a1b8e93
I am using DifferentialEquations version 4.5.0