# DiffEq skips over event in the first `solve` but works after `reinit!`

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.

``````using DifferentialEquations
using Plots
using BenchmarkTools

du = (sqrt(2) * 230 * sin(2*pi*50.0*t) - u * p) / 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 == 1
u
else
return 1
end

end

function affect_off!(integrator)
# println("applied affect off at t=\$(integrator.t)")
integrator.p = 1e10
integrator.p = 0
end

# test if diode needs to be turned on
function condition_on(u, t, integrator)
if integrator.p == 0
# @show integrator.t
# @show (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator) + u * 100)) - 1
return (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator) + u * 100)) - 1
else
return -1
end
end

function affect_on!(integrator)
# println("applied affect on at t=\$(integrator.t)")
integrator.p = 100
integrator.p = 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

``````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:

``````integrator.t = 0.03984942934492531
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator)) + u * 100)) - 1 = -14.524901788326165
integrator.t = 0.04421192005119001
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator)) + u * 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.

``````integrator = init(prob, solver,
callback=cb_set,
);
``````

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> maximum(sol.u)

1-element Vector{Float64}:
2.0098042506364657
``````

Please let me know if I missed something obvious. Thank you for your time.

I tested the code with `solver = Rosenbrock23()`, and it works fine upon the first `solve!`.

I’m using DifferentialEquations v7.6.0.

First thing, did you try lowering the tolerances?

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

`repeat_nudge=1//100` to something like `repeat_nudge=1//1000`. But anyways,

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

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.

Lowering the tolerance and using a different solver works. It’s just that the total time is greater than a fixed-stepping method.

How does the accuracy compare? It’s really easy to be faster if you are a lot less accurate.

And what about `repeat_nudge=1//10`?

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

``````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?

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?

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.

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.

``````solver = Rosenbrock23()

integrator = init(prob, solver,
callback=cb_set,
abstol=1e-3,
reltol = 1e-5,
);
``````

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.

Thanks again! You are right, `dt` keeps increasing. Setting `dtmax` solves the problem. Below is the solver and tolerance I’m using now:

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

``````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 (Code Optimization for Differential Equations · DifferentialEquations.jl) 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!

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.

1 Like

Noted. Thank you!!

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.

1 Like

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.

1 Like