# Local Error Estimation for DAEs
**URL:** https://discourse.julialang.org/t/local-error-estimation-for-daes/121159
**Category:** Modelling & Simulations
**Tags:** diffeq
**Created:** [October 10, 2024, 4:44pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159 "2024-10-10T16:44:11Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)
#### Post date: [October 10, 2024, 4:44pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/1 "2024-10-10T16:44:11Z")
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:
```Julia
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!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 10, 2024, 5:47pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/2 "2024-10-10T17:47:55Z")
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.
> [@aaliqureshi](#):
>
> - 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).
Because if not, then the solution is not guaranteed to be accurate. There are cases where it can be infinitely off.
> [@aaliqureshi](#):
>
> 1. Given that algebraic variables are solved exactly, what are the implications of disabling error checking for these variables on the global error?
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.
---
### Author:  [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)
#### Post date: [October 10, 2024, 8:45pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/3 "2024-10-10T20:45:27Z")
Hi Chris, thank you for the response. You mentioned that Rodas methods do not consider the error in algebraic variables.
> [@ChrisRackauckas](#):
>
> 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.
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?
> [@ChrisRackauckas](#):
>
> Because if not, then the solution is not guaranteed to be accurate. There are cases where it can be infinitely off.
If the numerical method is A-stable, wouldn’t the errors remain bounded?
Thanks!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 10, 2024, 9:25pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/4 "2024-10-10T21:25:37Z")
> [@aaliqureshi](#):
>
> 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?
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.
> [@aaliqureshi](#):
>
> If the numerical method is A-stable, wouldn’t the errors remain bounded?
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.
---
### Author:  [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)
#### Post date: [October 10, 2024, 9:40pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/5 "2024-10-10T21:40:37Z")
Thank you, Chris!
---
### Author:  [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)
#### Post date: [October 11, 2024, 2:04am UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/6 "2024-10-11T02:04:57Z")
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.
```julia
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,
> [@ChrisRackauckas](#):
>
> 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.
Does that mean the inclusion of errors from algebraic variables is an ‘extra step’ here to maintain numerical accuracy/stability?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 11, 2024, 2:47am UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/7 "2024-10-11T02:47:52Z")
> 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
---
### Author:  [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)
#### Post date: [October 11, 2024, 3:09pm UTC](https://discourse.julialang.org/t/local-error-estimation-for-daes/121159/8 "2024-10-11T15:09:15Z")
Thanks!