Simple DAE does not solve

Specific Domains Modelling & Simulations
,
1

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?

2

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.

3

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.