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

Specific Domains Numerics
, ,
1

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?

Thanks for your help,
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
2

I think it could help to use save_position

3

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.

4

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

5

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.

6

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.

7

@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.

8

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.

9

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.

10

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.

11

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.

12

Thanks for your answer.
My intuition was led by mostly dealing with continuous-time stochastic processes, which sample functions belong to a Skorokhod space, so I was just thinking about right-continuous processes. But of course, this might not be the best in every situation.

Thanks for pointing me to the additional continuity argument.