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

**URL:** <https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370>\
**Category:** Numerics\
**Tags:** question, sciml, differentialequation\
**Created:** [October 25, 2023, 8:50am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370 "2023-10-25T08:50:06Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 25, 2023, 8:50am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/1 "2023-10-25T08:50:06Z")

</div>

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:

```julia
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.

```julia
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`.

```julia
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 25, 2023, 10:31am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/2 "2023-10-25T10:31:19Z")

</div>

I think it could help to use `save_position`

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 25, 2023, 10:52am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/3 "2023-10-25T10:52:02Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 25, 2023, 11:22am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/4 "2023-10-25T11:22:28Z")

</div>

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](https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl) 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](https://arxiv.org/abs/1504.06873), arXiv preprint, 2015

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [October 25, 2023, 12:14pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/5 "2023-10-25T12:14:06Z")

</div>

> [@VincentW](#):
>
> ```julia
> # 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)
> 
> ```

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [October 25, 2023, 12:15pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/6 "2023-10-25T12:15:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 25, 2023, 4:06pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/7 "2023-10-25T16:06:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 30, 2023, 2:51pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/8 "2023-10-30T14:51:36Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [October 30, 2023, 3:10pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/9 "2023-10-30T15:10:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 30, 2023, 3:39pm UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/10 "2023-10-30T15:39:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 31, 2023, 8:32am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/11 "2023-10-31T08:32:31Z")

</div>

> [@VincentW](#):
>
> 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.

> **[Solution Handling · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#Interpolations-and-Calculating-Derivatives)**
>
> Documentation for DifferentialEquations.jl.

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.

---

<div class="post-metadata">

**Author:** ![VincentW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vincentw/32/202719_2.png) [@VincentW](https://discourse.julialang.org/u/VincentW)\
**Post date:** [October 31, 2023, 8:58am UTC](https://discourse.julialang.org/t/missing-precision-in-the-solution-of-combinations-of-ode-and-jump-processes-using-differentialequations-jl-and-jumpprocesses-jl/105370/12 "2023-10-31T08:58:48Z")

</div>

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.