Local Error Estimation for DAEs

Specific Domains Modelling & Simulations
1

I’m inquiring about the error estimation approach used for solving Differential Algebraic Equations (DAEs) in Julia.

The current approach considers all variables (state and algebraic) in the estimation of local error. However, this seems to impose (unnecessary) smoothness restrictions on algebraic variables.

For example, when solving DAEs with discontinuities without re-initialization, the solver reports instability at the first step following the discontinuity. This instability arises from inconsistent initial conditions, which causes large jumps in some variables. I observed that relaxing the tolerance for algebraic variables allowed the solver to proceed (see MWE).

My questions are:

  1. Why are algebraic variables involved in local error estimation during numerical integration of DAEs? Since numerical method isn’t directly applied to algebraic equations, why are their solutions used in estimating Local Truncation Error (LTE).
  2. Given that algebraic variables are solved exactly, what are the implications of disabling error checking for these variables on the global error?

MWE:

using OrdinaryDiffEq

function fx!(du,u,p,t)
    x,y,z = u

    du[1] = y
    du[2] = -x
    du[3] = x^2 + y^2 - z^2 -1
end


# inconsistent initial condition
u0 = [2,0,1]

# inconsistent initial conditions but close to consistent initial condition
# u0 = [2, 0, 1.7]


# this is consistent initial condition
# u0 = y1 = [1.999999995, -9.999999975e-5, 1.7320508046821261]

p = ()
M = zeros(3,3)
M[1,1] = M[2,2] = 1

dt = 50e-6
t_span = (0,1)

prob0 = ODEFunction(fx!, mass_matrix = M)
prob = ODEProblem(prob0, u0,t_span, p)

#define tolerance
tol = [1e-7, 1e-7, 1]

# solve with uniform reltol for all variables (this fails)
sol_1 = solve(prob, Rodas3(), adaptive = true, reltol = 1e-7, initializealg = NoInit())

# solve with relaxed tolerance for algebraic variables (this works)
sol_2 = solve(prob, Rodas3(), adaptive = true, reltol = tol, initializealg = NoInit())

Thank you!

2

Note that it’s dependent by method. On the Rodas methods, the algebraic condition error is not used in the error estimation. Rosenbrock23 directly enforces it. It uses an embedded estimate which compares state vectors. Rodas5Pr is somewhat of a hybrid via the interpolation.

You could however pass a vector of tolerances and have a very high tolerance on the algebraic equations.

Because if not, then the solution is not guaranteed to be accurate. There are cases where it can be infinitely off.

If you’re using Rodas5P, it should do exactly that. Which is the reason for Rodas5Pr for more robustness.

The implication is that you can have numerical drift due to floating point truncation error, and that can cause the initialization to not be satisfied in the next step, and that can then have a compounding effect.

1 Like
3

Hi Chris, thank you for the response. You mentioned that Rodas methods do not consider the error in algebraic variables.

But in the MWE presented above, I used Rodas3 solver and it fails when a tighter tolerance is used for algebraic variables and works fine when the said tolerance is increased.

Also, I am not sure I understood what you meant here. Could you please elaborate?

If the numerical method is A-stable, wouldn’t the errors remain bounded?

Thanks!

4

The method “by definition” will solve the algebraic variables “exactly”, so then the error is not considered in the error estimate itself. But if it cannot project it essentially is a solver error.

No, A-stable is a property just a linera property (and on ODEs). Stiffly stable and B-stability should be considered on DAEs. Rosenbrock23 for example can technically solve DAEs in mass matrix form without a change, but it will have no error control on the algebraic conditions without changing the error estimator.

1 Like
5

Thank you, Chris!

6

Hi Chris, I looked into the implementation of Rosenbrock23 solver and found that the errors in algebraic variables are indeed included in the error estimate.

        if mass_matrix !== I
            algvar = reshape(cache.algebraic_vars, size(u))
            @.. broadcast=false atmp=ifelse(algvar, fsallast, false) /
                                     integrator.opts.abstol
            integrator.EEst += integrator.opts.internalnorm(atmp, t)
        end

As you said before,

Does that mean the inclusion of errors from algebraic variables is an ‘extra step’ here to maintain numerical accuracy/stability?

7

Does that mean the inclusion of errors from algebraic variables is an ‘extra step’ here to maintain numerical accuracy

yes. for singular mass matrices, you probably should prefer Rodas23W over Rosenbrock23

2 Likes
8

Thanks!