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

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

```
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> maximum(sol.u)
1-element Vector{Float64}:
2.0098042506364657
```

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