# DiffEq skips over event in the first \`solve\` but works after \`reinit!\`
**URL:** https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935
**Category:** Modelling & Simulations
**Created:** [November 28, 2022, 5:08am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935 "2022-11-28T05:08:03Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [November 28, 2022, 5:08am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/1 "2022-11-28T05:08:03Z")
The code below simulates a series resistor and inductor supplied by an AC source through a diode. The diode turns off when the current is \<=0 and turns back on when the voltage across the diode is \>=1.
# Skipping over events
My first issue with the code is that it skips over a turn-on event the first time it’s run. The second run yields the expected result.
```julia
using DifferentialEquations
using Plots
using BenchmarkTools
function I_load!(du, u, p, t)
du[1] = (sqrt(2) * 230 * sin(2*pi*50.0*t) - u[1] * p[1]) / 0.5
end
u0 = [0.0]
p = [100, 1] # resistance, diode flag
tspan = (0.0, 60.0E-3)
# test if diode needs to be turned off
function condition_off(u, t, integrator)
# println("test affect off at $(integrator.t), u=$(integrator.u)")
if integrator.p[2] == 1
u[1]
else
return 1
end
end
function affect_off!(integrator)
# println("applied affect off at t=$(integrator.t)")
integrator.p[1] = 1e10
integrator.p[2] = 0
end
# test if diode needs to be turned on
function condition_on(u, t, integrator)
if integrator.p[2] == 0
# @show integrator.t
# @show (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
return (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
else
return -1
end
end
function affect_on!(integrator)
# println("applied affect on at t=$(integrator.t)")
integrator.p[1] = 100
integrator.p[2] = 1
end
cb_off = ContinuousCallback(condition_off, nothing, affect_off!) # only when condition goes negative
cb_on = ContinuousCallback(condition_on, affect_on!, nothing) # only when condition goes positive
cb_set = CallbackSet(cb_off, cb_on);
prob = ODEProblem(I_load!, u0, tspan, p)
solver = AutoTsit5(Rosenbrock23())
integrator = init(prob, solver,
callback=cb_set,
);
sol = solve!(integrator)
plot(sol)
```
Doing a `reinit!` yields the expected solution
```julia
reinit!(integrator, u0)
prob.p[:] = [100, 1.0]; # restore parameters
sol = solve!(integrator)
plot(sol.t, sol')
```
If you uncomment the two lines in function `condition_on`, the following output can be seen:
```julia-auto
integrator.t = 0.03984942934492531
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator))[1] + u[1] * 100)) - 1 = -14.524901788326165
integrator.t = 0.04421192005119001
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator))[1] + u[1] * 100)) - 1 = 51.91978477325262
```
This means that the function reports a rising zero crossing, but the solve did not enter the time interpolation to detect the zero crossing. I spent several hours in troubleshooting my event functions but didn’t see obvious issues here.
# Trapezoid() gives wrong results - Stiffness?
Another issue about the solver. If the solver is changed to `Trapezoid()` or `AutoTsit5(Trapezoid())`, the results are wrong.
`AutoTsit5(Trapezoid())` plot using `plot(sol.t, sol')`
I speculate that the Trapezoid method also skips over events. With a really tiny step size, the solution is expected.
```julia
integrator = init(prob, solver,
callback=cb_set,
dt=5e-7, adaptive=false
);
```
This step size is obviously too small to be practical. The same problem is solved by a hand-baked trapezoid method in a textbook using a step size of 5e-5 (although that example is used to demonstrate numerical chattering). But if I use `dtmax=5e-5, adaptive=true`, the last turn-on event is again skipped, producing the same result as in the first plot.
# Different Plots
Also, `plot(sol.t, sol')` and `plot(sol)` gives seemingly unrelated plots.
`AutoTsit5(Trapezoid())` plot using `plot(sol)`:
```julia
julia> maximum(sol.u)
1-element Vector{Float64}:
2.0098042506364657
```
Please let me know if I missed something obvious. Thank you for your time.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [November 28, 2022, 7:48pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/2 "2022-11-28T19:48:09Z")
I tested the code with `solver = Rosenbrock23()`, and it works fine upon the first `solve!`.
I’m using DifferentialEquations v7.6.0.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 29, 2022, 7:22pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/3 "2022-11-29T19:22:22Z")
First thing, did you try lowering the tolerances?
> [@cuihantao](#):
>
> `integrator.p[1] = 1e10`
That’s going to be ill-conditioned, so you’ll need a lower tolerance solution in order to resolve it well. My guess is that a lower tolerance is all you need, and that at higher tolerances you may need to adjust
[https://diffeq.sciml.ai/stable/features/callback\_functions/#SciMLBase.ContinuousCallback](https://diffeq.sciml.ai/stable/features/callback_functions/#SciMLBase.ContinuousCallback)
`repeat_nudge=1//100` to something like `repeat_nudge=1//1000`. But anyways,
```julia
integrator = init(prob, solver,
callback=cb_set,
abstol = 1e-8, reltol = 1e-8
);
```
is probably all that’s needed. And thn I wouldn’t use `Rosenbrock23` for those kinds of tolerances.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [November 29, 2022, 7:24pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/4 "2022-11-29T19:24:13Z")
Thanks a lot, Chris. I was aware of the ill condition but did not understand that a low tolerance solution is achievable in such a condition. I will test it and report back.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 2, 2022, 6:17pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/5 "2022-12-02T18:17:15Z")
Lowering the tolerance and using a different solver works. It’s just that the total time is greater than a fixed-stepping method.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 2, 2022, 6:40pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/6 "2022-12-02T18:40:04Z")
How does the accuracy compare? It’s really easy to be faster if you are a lot less accurate.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 2, 2022, 7:21pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/7 "2022-12-02T19:21:11Z")
And what about `repeat_nudge=1//10`?
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 2, 2022, 8:13pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/8 "2022-12-02T20:13:28Z")
It does cut the total time by a lot. From `destats`, I can tell the number of root-finding calls has been reduced. Currently, I’m using
```julia
solver = Rodas4P()
integrator = init(prob, solver,
callback=cb_set,
abstol=1e-3,
reltol = 1e-5,
);
```
with `repeat_nudge=1//10`.
The total time is now 116.616 μs (643 allocations: 30.81 KiB). The custom integrator takes about 2 us (without allocation). That’s not an apples-to-apples comparison but still two orders apart. I couldn’t identify allocations in my code. How could I further reduce allocations?
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 2, 2022, 8:16pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/9 "2022-12-02T20:16:35Z")
The DiffEq results are more accurate right after the event, but I’m willing to trade reasonable accuracy for speed.
Below shows the results using a custom integrator (based on the trapezoidal rule) and using Rodas4P. What settings can I change to further improve the speed?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 2, 2022, 8:26pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/10 "2022-12-02T20:26:13Z")
If you want to trade accuracy for speed, raise your tolerances. `Rodas4P` is also usually not the best solver. `Rosenbrock23` will be probably be better with loose tolerances, and `Rodas5P` will usually be better at stricter tolerances.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 2, 2022, 8:30pm UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/11 "2022-12-02T20:30:31Z")
I understand that a lower tolerance leads to a higher number of iterations and thus slower speed.
Using `Rosenbrock23`, the results are not correct. It skips over the events. That’s why I posted the original problem.
`repeat_nudge=1//10` and the following is what I used to get the figure above.
```julia
solver = Rosenbrock23()
integrator = init(prob, solver,
callback=cb_set,
abstol=1e-3,
reltol = 1e-5,
);
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2022, 1:09am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/12 "2022-12-03T01:09:23Z")
If it’s missing an event, you may need to decrease it. Can you have it print out `integrator.dt` in the `affect!`? I think this may be a case where you have a bit of an adverse effect because the ODE might be simple enough that the error is essentially zero for any step, so the `dt` always grows, and so it grows so large that numerical errors become easy. In which case, fixing `maxdt` needs to be done. Show `integrator.dt` to see if that’s the case.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 3, 2022, 3:28am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/13 "2022-12-03T03:28:33Z")
> [@ChrisRackauckas](#):
>
> Can you have it print out `integrator.dt` in the `affect!`?
Thanks again! You are right, `dt` keeps increasing. Setting `dtmax` solves the problem. Below is the solver and tolerance I’m using now:
```julia
solver = Rodas5P()
integrator = init(prob, solver,
callback=cb_set,
abstol=1e-4,
reltol = 1e-5,
dtmax=1/60,
);
```
I need to build a larger problem to compare DiffEq’s solver time over a manual fixed-step variant of trapezoidal rule. Currently, DiffEq’s time is 122 us and the manual solver just takes 4 us. It’s probably because the manually written solver is non-allocating, and variables happen to fit into the cache.
The current `destats` are
```julia
DiffEqBase.DEStats
Number of function 1 evaluations: 384
Number of function 2 evaluations: 0
Number of W matrix evaluations: 39
Number of linear solves: 312
Number of Jacobians created: 35
Number of nonlinear solver iterations: 0
Number of nonlinear solver convergence failures: 0
Number of rootfind condition calls: 843
Number of accepted steps: 35
Number of rejected steps: 4
```
There is a total of 645 allocations: 32.42 KiB. The optimization example here ([https://diffeq.sciml.ai/stable/tutorials/faster\_ode\_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation](https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation)) only has 41 allocations. Is my allocation number normal?
Also, since my actual problem is going to be large, I won’t be able to use the optimization techniques for small matrices.
Thank you again!
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 3, 2022, 3:33am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/14 "2022-12-03T03:33:07Z")
You can use `saveat` to get preallocated output storage. That said, if your actual problem will be large, benchmarking small problems isn’t that useful. The larger problem will probably be stiffer, and require more of the machinery that the fixed timestep version is missing.
---
### Author:  [@cuihantao](https://discourse.julialang.org/u/cuihantao)
#### Post date: [December 3, 2022, 3:37am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/15 "2022-12-03T03:37:19Z")
Noted. Thank you!!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2022, 3:41am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/16 "2022-12-03T03:41:27Z")
> [@cuihantao](#):
>
> There is a total of 645 allocations: 32.42 KiB. The optimization example here ([Code Optimization for Differential Equations · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation)) only has 41 allocations. Is my allocation number normal?
A lot of those allocations are just saving the time series. I think there’s some more in the callback code that I need to eliminate though.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 3, 2022, 3:42am UTC](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/17 "2022-12-03T03:42:06Z")
Specifically, for large systems all the time is in jacobian solves and function evaluation. Allocating a length 3 array is roughly 1/3d as expensive as allocating a length 100 array, but the jacobians will be thousands of times slower.