I tried the very most basic forced DAE I could think of and implemented an RC-circuit:
using Sundials
using OrdinaryDiffEq
using Test
u(t) = t < 0.5 ? 0.0 : 1.0
p = (R=1.0, C=1.0)
righthandside = (out, dx, x, p, t) -> begin # simple RC circuit DAE
# x: [v_C, v_R, i_C, i_R]
out[1] = 1 / p.C * x[3] - dx[1] # dv_C/dt = i_C/C
out[2] = x[2] - p.R * x[4] # v_R - R*i_R = 0 (Standard Ohm's law)
out[3] = x[1] + x[2] - u(t) # v_C + v_R - u(t) = 0
out[4] = x[3] - x[4] # i_C - i_R = 0 (Series current equality)
return out
end
prob = DAEProblem(
righthandside,
zeros(4), # initial state derivative
zeros(4), # initial state
(0.0, 10.0),
p;
differential_vars=[true, false, false, false],
)
integrator = init(prob, IDA(); reltol=1e-6, abstol=1e-6, initializealg=OrdinaryDiffEq.BrownFullBasicInit())
for t in range(0.0; stop=10.0, step=0.1)
step!(integrator, 0.1, true)
end
@test OrdinaryDiffEq.SciMLBase.successful_retcode(sol)
@test sol.t[end] ≈ 10.0
However, I get an error from the solver.
julia> include("dae.jl")
[ERROR][rank 0][/workspace/srcdir/sundials/src/idas/idas.c:5732][IDAHandleFailure] At t = 0.5 and h = 1.06658388718261e-17, the error test failed repeatedly or with |h| = hmin.
Test Passed
Is stepping the integrator manually to the discontinuity and using a callback not the same?
No because it needs to reinit after that. u(t) having a first order discontinuity means that your solution is discontinuous. For ODEs, having a discontinuity like that in the right hand side means your derivative is discontinuous, but here you have a DAE with a discontinuity in your algebraic equations which means that x[2](t) is actually discontinuous. This can make the time stepping not converge because stepping near that point means that delta x[2](t) / dt does not go to zero as dt->0. This means the error control does not converge to zero. What needs to happen then is a discontinuity needs to be defined at that point and the solver needs to reinit! after the discontinuity, which is done automatically in the callback handling.
Or using ModelingToolkit, x[2] is eliminated from the solving process which stabilizes the solving. I highly recommend MTK for any DAE because it just handles a lot of these things: lots of DAEs are not necessarily numerically solvable in the form people first write them.
As a follow up. If an MTK model cannot be converted to an ODE, I still have this problem if there is a discontinuity involved? Which I always have in a hybrid system, i.e. if u(t) would be sampled.
The docs state that this exact problem “changing an algebraic variable” might not even work with callbacks either: Event Handling and Callback Functions · DifferentialEquations.jl
So the solution: Use MTK. 
MTK’s solution is to ensure that the callback has an implicit definition that forces satisfaction. The way that this is done is to include the algebraic equations as part of an implicit nonlinear system for the callback. So let’s say your callback is:
function affect!(integrator)
intregrator.u[1] += 1
end
For an ODE, that’s fine. But for a DAE, you also need to ensure the algebraic conditions are satisfied, which if you change x[1] then the 3rd condition here is not satisfied because x[2] needs to update. So what the MTK symbolic callbacks would do instead is you say:
x1 ~ Pre(x1) + 1
Which it then automatically expands your case here to:
x1 ~ Pre(x1) + 1
0 ~ x2 - p.R * x4
0 ~ x1 + x2 - u(t)
0 ~ x3 - x4
it would then simplify this, x4 ~ x3, x3 ~ x2 / p.R, x2 ~ u(t) - x1, or in other words “what you truly meant” by intregrator.u[1] += 1 is actually:
function affect!(integrator)
integrator.u[1] += 1
integrator.u[2] = u(integrator.t) - integrator.u[1]
integrator.u[3] = integrator.u[2] / p.R
integrator.u[4] = integrator.u[3]
end
which is then an explicit update system that satisfies the DAE conditions and requires no separate reset.
Then there’s tools like iflifting that then take into account the discontinuity, so it can split Pre(x1) is pre-discontnuity and x1 is post, meaning that by splitting this you can get solvable equations.