# Missing precision in the solution of combinations of ODE and jump processes using DifferentialEquations.jl and JumpProcesses.jl

Hey everyone,

For a toy model, I tried to simulate an ODE together with a jump process, which jump rate at a certain time is dependent on the ODE at that time. (Motivated by the growth of a tumor and a death process, where the probability to die increases by tumor size.)

Using DifferentialEquations.jl and JumpProcesses.jl this can be done very cleanly as shown in this MWE:

``````using DifferentialEquations
using JumpProcesses
using Plots
using Random

# set seed for reproducibility
Random.seed!(123)

# define the drift function
function ODE!(du, u, p, t)
if (u[2] == 0)
du[1] = p.Î˛*u[1]
else
du[1] = 0
end
end

# define the jump-rate functions
jump_rate(u, p, t) = p.Î»*u[1]^2

# define the jump functions
function jump_affect!(integrator)
if (integrator.u[2] == 0.0) # still alive
integrator.u[2] += 1
end
nothing
end

# define initial conditions and parameters
u0 = [0.05, 0.0]
p = (Î˛ = 0.3, Î» = 0.1,)
tspan = (0.0, 30.0)

# define the ODE problem
prob = ODEProblem(ODE!, u0, tspan, p)

# define the jump problem
vrj = VariableRateJump(jump_rate, jump_affect!)
prob = JumpProblem(prob, Direct(), vrj)

sol = solve(prob, Tsit5());
``````

However, when I access the solution, I run into imprecisions and inconsistent results.

I wanted to find out the death event time and assumed the solver stores the times, where the process jumps. So I searched for them via the array interface of the solution object.

``````jump_index = findfirst([tuples(sol)[i][1][2] == 1.0 for i in 1:length(tuples(sol))]);
jump_time = tuples(sol)[jump_index][2];
``````

But when I call the solution object as a function of that time, I get a different result than the one that is stored in the `sol.u` array for that timepoint. I.e. the second line of the following code block will evaluate to `false`.

``````sol.t[jump_index] == jump_time

sol.u[jump_index] == sol(jump_time)
``````

I am aware that using the solution object as a function of a timepoints, just returns an approximated/interpolated value of the solution.
However, it seems inconsistent to me that this also happens for the stored values. In my opinion they should agree with each other at those timepoints.

The explanation I could imagine to make sense here, is that the solution object stores the timepoints not to full precision, but cuts them after a certain number of decimal places.

Still, I am not sure whether this is the reason and find it nevertheless inconsistent and not ideal especially when working with jump processes, where I am interested in the exact jump times of the simulation.

Does any one know why this behavior occurs?
Is there any reason why this behaviour canâ€™t be circumvented or is even desirable?

Vincent

Ps: I hope the above MWE is understandable. If not please let me know how to improve on it or the title of the post.

1 Like

I think it could help to use `save_position`

Do you mean for saving the jump times?
I guess this could then be achieved, however this behavior still seems to be inconsistent, for me.

I think your system has an analytical solution. Anyway I treated more or less the same in GitHub - rveltz/PiecewiseDeterministicMarkovProcesses.jl: Piecewise Deterministic Markov Processes in Julia which uses a different algo from DE.

It is based on R. Veltz, A new twist for the simulation of hybrid systems using the true jump method, arXiv preprint, 2015

If you all you want is the first time that the jump occurs you can use `terminate!` to end the simulation then

``````function jump_affect!(integrator)
if (integrator.u[2] == 0.0) # still alive
integrator.u[2] += 1
end
terminate!(integrator)
nothing
end
``````

Then `sol.t[end]` should be this time.

I agree that something isnâ€™t quite right when evaluating the solution object via interpolation, and I suspect it has to do with the ODE interpolation being used across the jump time. Hopefully @ChrisRackauckas can help with that. If you use the modified `jump_affect!` above the interpolation `sol(sol.t[end])` still seems to have an issue.

So as @rveltz suggested, I would suggest manually specifying times you want to save the solution at, and using `sol.u` to get the values, for now.

@rveltz @isaacsas Thanks for your suggestion how to get death time and corresponding simulation values. I am doing it for now using `terminate!`.

Still, I hope that someone might fix the issue with the discrepancy between interpolation solution and stored values.

Using `terminate!` I encountered that this issue might be even worse than just simple imprecisions at the jump time.

If one uses `terminate!` as suggested by @isaacsas

``````function jump_affect!(integrator)
if (integrator.u[2] == 0.0) # still alive
integrator.u[2] += 1
end
terminate!(integrator)
nothing
end
``````

Then the solution object will store the solution at jump-time twice, once with the jump and once without and then the process ends.

As pointed out before, if one then calls `sol(jumptime) it will return the state of the process, without the jump.

Additionally, for any time `t>jumptime`, the interpolation `sol(t)` will return the same value of the process as before the jump happened. So the jump process will have the value 0.0, although a jump occured and was stored in the solution object.

it makes sense that the process stays constant since the integrator terminated after the jump.
However, it does not seem right, that the jump process stays constant at 0.0 and not at 1.0.

You should never evaluate a solution object beyond the final time. The interpolations are only meant to be used within the time interval of solving. If you need solutions beyond that time donâ€™t use terminate.

I agree that one should not evaluate a solution object beyond the final time. I was more generally playing around with the solution object to understand its behavior better.

Still, I find this behavior somehow unintuitive and think it would be better/more precise if the value after the jump would be returned.

Intuitive depends on who youâ€™re talking to. Itâ€™s simply undefined because thereâ€™s two values of the solution at that `t`. Undefined canâ€™t have a â€śbetter/more preciseâ€ť answer since both are true. Given that itâ€™s undefined, thereâ€™s an argument `continuity` for controlling which solution you get.

You can use `sol(t, continuity = :right)` if you want the right-continuous value. There are just as many cases (well, arguably more) where this choice is less intuitive which is why the default is set as it is.

Thanks for pointing me to the additional `continuity` argument.