Lecture 2 – Simulation of differential-algebraic equations

Simulation of DAE models and index reduction

Erik Frisk
erik.frisk@liu.se

Department of Electrical Engineering

Linköping University

May 21, 2024

1 / 50
Simple circuit model, index 1
iL
e1 : u0 = f (t)
u1
e2 : u1 = R1 i1
L uL

R1
i0 e3 : u2 = R2 i2
u2
i1 duc
e4 : iC = C
i2
R2 dt
diL
e5 : uL = L
U0 dt
iC
e6 : i0 = i1 + iL
uC C
e7 : i1 = i2 + iC
e8 : u0 = u1 + uC
e9 : uL = u1 + u2
e10 : uC = u2
x1 = (uc , iL ), x2 = (u2 , i2 , u0 , u1 , uL , i1 , iC , i0 )

2 / 50
Computational form of model
duc 1
e4 : = iC
dt C
e1 : u0 = f (t) diL 1
e2 : u1 = R1 i1 e5 : = uL
dt L
e3 : u2 = R2 i2
duc e10 : u2 := uC
e4 : iC = C
dt 1
e3 : i2 := u2
diL ⇒ R2
e5 : uL = L
dt e1 : u0 := f (t)
e6 : i0 = i1 + iL e8 : u1 := u0 − uC
e7 : i1 = i2 + iC e9 : uL := u1 + u2
e8 : u0 = u1 + uC 1
e9 : uL = u1 + u2 e2 : i1 := u1
R1
e10 : uC = u2 e7 : iC := i1 − i2
e6 : i0 := i1 + iL
3 / 50
Simple circuit model, index > 1
iC
e1 : u0 = f (t)
e2 : u1 = R1 i1
u1 C uC

R1
i0 e3 : u2 = R2 i2
u2
i1 duc
e4 : iC = C
i2
R2 dt
diL
U0
e5 : uL = L
iL dt
e6 : i0 = i1 + iC
uL L e7 : i1 = i2 + iL
e8 : u0 = u1 + uL
e9 : uC = u1 + u2
e10 : uL = u2
x1 = (uC , iL ), x2 = (u2 , i2 , u0 , u1 , uL , i1 , iC , i0 )

4 / 50
Differential index

F (t, y , ẏ ) = 0

Definition
The minimum number of times the DAE has to be differentiated with
respect to t to be able to determine ẏ as a function of t och y is called
the (differential-) index of the DAE.

index might be solution dependent, uniform index

There are several types of index, the above is called differential index.
Perturbation index
variants of the above (see paper)
Anyhow: index is a measure how far from an ODE the DAE is.

5 / 50
Outline

Simulation of high index DAE:s, key problems

Simulation of index 1 DAE:s

State-space method
ϵ-embedding
BDF

Index reduction
Index reduction by differentiation
Drift stabilization

6 / 50
Simulation of DAE with index > 1

Aẏ (t) + By (t) = g (t), linear, constant coefficients

A(t)ẏ (t) + B(t)y (t) = g (t), linear, time-varying coefficients
F (ẏ , y , t) = 0, general DAE

7 / 50
Simulation of DAE with index > 1

Aẏ (t) + By (t) = g (t), linear, constant coefficients

A(t)ẏ (t) + B(t)y (t) = g (t), linear, time-varying coefficients
F (ẏ , y , t) = 0, general DAE

Start with a linear DAE with constant coefficients, this is sufficient to

illustrate the main reasons why high index problems are difficult
Consider the index-3 problem

x(t) = g (t)
ẋ = y
ẏ = z

which has the solution x(t) = g (t), y (t) = ġ (t), and z(t) = g̈ (t).
7 / 50
Backward Euler, fix step length
A backward Euler on the problem gives the equations

xn = gn
xn − xn−1
yn =
h
yn − yn−1
zn =
h
Direct substitutions of xn and yn give

xn = gn
gn − gn−1
yn =
h
gn − 2gn−1 + gn−2
zn =
h2

8 / 50
Backward Euler, fix step length
A backward Euler on the problem gives the equations

xn = gn
xn − xn−1
yn =
h
yn − yn−1
zn =
h
Direct substitutions of xn and yn give

xn = gn
gn − gn−1
yn = = ġn + O(h)
h
gn − 2gn−1 + gn−2
zn = = g̈n + O(h)
h2

8 / 50
Backward Euler, fix step length
A backward Euler on the problem gives the equations

xn = gn
xn − xn−1
yn =
h
yn − yn−1
zn =
h
Direct substitutions of xn and yn give

xn = gn
gn − gn−1
yn = = ġn + O(h)
h
gn − 2gn−1 + gn−2
zn = = g̈n + O(h)
h2
This looks great!

x(tn ) − xn = 0, y (tn ) − yn = O(h), z(tn ) − zn = O(h)

8 / 50
Fix step length

Theorem
If a k step BDF (k < 7) is applied to

Aẏ (t) + By (t) = g (t)

the solution will be O(hk ) after max (m − 1)k + 1 steps.

This would indicate that an algorithm with fixed step length would work
fine also for high index.
Unfortunately, this breaks down for variable step length integrators.

9 / 50
Constant DAE, variable step length
Now assume variable step length, i.e., the algorithm becomes

xn = gn
xn − xn−1
yn =
hn
yn − yn−1
zn =
hn

10 / 50
Constant DAE, variable step length
Now assume variable step length, i.e., the algorithm becomes

xn = gn
xn − xn−1
yn =
hn
yn − yn−1
zn =
hn
Now, examining the error in zn we obtain
 
1 hn−1
zn − gn′′ = · · · = − 1 g̈n + O(h)
2 hn

10 / 50
Constant DAE, variable step length
Now assume variable step length, i.e., the algorithm becomes

xn = gn
xn − xn−1
yn =
hn
yn − yn−1
zn =
hn
Now, examining the error in zn we obtain
 
1 hn−1
zn − gn′′ = · · · = − 1 g̈n + O(h)
2 hn

This means that the error will diverge(!) with decreasing hn , O(hn−1 )
(index > 2). The error in y will be O(1) (index 2).
One of the exercises is to show . . . in the expression.
Hint: Taylor expansion around t = tn .
10 / 50
Constant DAE, variable step length, cont.

One result for variable step length

Theorem
If a k step BDF (k < 7) is applied on

Aẏ (t) + By (t) = g (t)

and the ratio between successive step lengths are bounded, then the
q
solution will be O(hmax ) where q = min(k, k − m + 2).

This gives that an index 6 problem could be solved by a 6-step BDF,

or?
Which precondition is here questionable, and why?
One step backwards Euler with fixed step length is a recommended
approach for linear high-index problems with constant coefficients.

11 / 50
Linear, non-constant DAE:s
For a DAE
A(t)ẏ (t) + B(t)y (t) = g (t)
you can define local index at time t and global index via a transformation
to a canonical form by y (t) = H(t)z(t) and G (t) to

G (t)A(t)(H ′ (t)z(t) + H(t)ż(t)) + G (t)B(t)H(t)z(t) = G (t)g (t)

such that
   
I 0 ′ C (t) 0
G (t)A(t)H(t) = , G (t)A(t)H (t) + G (t)B(t)H(t) =
0 E 0 I

where E is nilpotent, i.e., the same canonical form as previously shown for
linear, constant, DAE:s

ż1 + C (t)z1 = G1 (t)g (t)

E ż2 + z2 = G2 (t)g (t)

12 / 50
Linear, non-constant DAE:s, cont.

If we can transform the DAE, all is well. Then the time variable
problem is no more difficult than the time invariant.
Problem: Finding the transformation matrices is not easy
What happens if you ”go with it” anyway with a fixed step length
BDF for a time variable system?
From Gear, Petzold:
If the local index is two, we may have a stability problem depending
on how fast the matrices change with time. If the local index is
larger than 2, we almost always have a stability problem.
Important to note: Stability problem, not an accuracy problem

13 / 50
What is the origin of the stability problem?
A DAE with local index m can be transformed into

(Pn EQn )z ′ + (Pn Qn )z = Pn qn

where Pn and Qn are time variable transformation matrices and the

constant matrix E is in the form (m = 3)
 
0 0 0
E = 1 0 0 , E m = 0
0 1 0

Apply BDF on the DAE

(Pn EQn )(zn − zn−1 ) + h(Pn Qn )zn = hPn qn

which then can be written as

(Pn EQn + hPn Qn )zn = Pn EQn zn−1 + hPn qn

14 / 50
Simplify to
(E + hI )Qn zn = EQn zn−1 + hqn
Now, solve for zn according to

zn = Qn−1 (E + hI )−1 EQn zn−1 + hQn−1 (E + hI )−1 qn = Sn zn−1 + un

The real solution satisfies

h2
z(tn ) = Sn z(tn−1 ) + un − Sn zn′′ (ξ)
2
With en = zn − z(tn ) we get the recursion

h2
en = Sn en−1 + Sn zn′′ (ξ)
2
Now we have a (recursive) expression for the simulation error, time to
analyze!

15 / 50
Expansion of the recursion gives

h2
en = Sn en−1 + Sn zn′′ (ξ)
2
where we have the explicit expression for the fault
 
2 n n n
h X Y
 Sj  zi′′ +
Y
en = Sj e0
2
i=1 j=i j=0

Three problems
z ′′ not bounded
e0 ̸= 0
Sj not bounded

16 / 50
Now we can start answering the questions, why does it work for constant
linear DAE:s (no time variable transformations):

Sn = S = Q −1 (E + hI )−1 EQ

Then we have that

S k = Q −1 ((E + hI )−1 E )k Q

Some basic algebra gives

 
0
 h−1 
 

−h−2 .. 
. 
−1
(E + hI ) E = 
 
 h−3 .. 
 . 

 .. 
 . 
h−3 h−2 h−1 0

which gives that S m = 0, i.e., S also is nilpotent of order m, same as E .

17 / 50
Back to the expression
 
n n n
h2 X Y
 Sj  zi′′ +
Y
en = Sj e0
2
i=1 j=i j=0

Here it is clear that for n > m, the second term disappears (that is not
influenced by the step length h)
After some thought you can derive the expression
m−2
h2 X i+1 ′′
en = S zn−i
2
i=0

Here we see why the limit for one step BDF is at index 1 problems (S
factors then contains at most h−1 )

18 / 50
Then why did it not work for time variable systems?
 
2 n n n
h X Y
 Sj  zi′′ +
Y
en = Sj e0
2
i=1 j=i j=0

Three problems
z ′′ not bounded
e0 ̸= 0
Sj not bounded
S2 S1 = Q2−1 (E + hI )−1 E Q2 Q1−1 (E + hI )−1 EQ1
| {z }
̸=I

For time variable/non-linear systems, the matrices are changing all the
time and ji=1 Si ̸= 0 for all j.
Q

19 / 50
Some conclusions

In the general case there is no methods (to my knowledge) for high index
problems in the form
A(t)ẏ + B(t)y = g
and even less for
F (t, ẏ , y ) = 0
Though, sometimes, it might work with a BDF with fixed step length.
Note that we have a stability issue, not accuracy (although this could
also happen, more about that next time)
Thus, it is not a solution to increase the order of the method or
taking shorter steps.
On the contrary, shorter steps might even make the situation worse
Classes of DAE:s

20 / 50
Semi-explicit DAE:s and a rule of thumb

If If
F (t, y , ẏ ) = 0
ẋ = f (x, y )
has index ν then
0 = g (x, y )
ẏ = u
has index ν then
0 = F (t, y , u)
ẋ = f (x, u̇)
index ν + 1.
0 = g (x, u̇)

index ν − 1.
Rule of thumb: The semi-explicit case behaves as the general but with a
higher index (and vice versa)

21 / 50
Outline

Simulation of high index DAE:s, key problems

Simulation of index 1 DAE:s

State-space method
ϵ-embedding
BDF

Index reduction
Index reduction by differentiation
Drift stabilization

22 / 50
Simulation of semi-explicit index 1 DAE:s

State-space method
ϵ-embedding
BDF (DASSL with variants) is a commonly used DAE solver. Can be
downloaded online.
Is described in detail in the nice book ”Numerical Solution of
Initial-Value Problems in Differential-Algebraic Equations” by K.E.
Brenan, S.L. Campbell and L.R. Petzold. The book is in our library
and I have uploaded the chapter on DASSL on the course web page
for the interested.
I will show the basic principles of DASSL at the end of this lecture.
I can also highly recommend the method descriptions of the
SUNDIALs solvers (https://sundials.readthedocs.io/).

23 / 50
Numerical code

Matlab have several solvers for DAE (ode15s, ode15i, ode23t)

Python has good ODE support in SciPy but no DAE solvers, I use the
package ODES https://scikits-odes.readthedocs.io/ which is
a wrapper around . . .
SUNDIALS (more on the next slide)
Julia – probably the environment with the most support for numerical
integration (but I think less support for DAE) with the package
DifferentialEquations.jl
(https://github.com/SciML/DifferentialEquations.jl)

24 / 50
SUNDIALS
SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers
Software library consisting of 6 different solvers, written in C.
https://computing.llnl.gov/casc/sundials

CVODE(S)
Solves IVP for ordinary differential equation (ODE) systems. Includes
sensitivity analysis capabilities (forward and adjoint).
IDA(S)
Solves IVP for differential-algebraic equation (DAE) systems. Includes
sensitivity analysis capabilities (forward and adjoint).
ARKode
Solves IVP ODE problems with additive Runge-Kutta methods,
including support for IMEX methods.
KINSOL
solves nonlinear algebraic systems.
Lots of functionality is not available in vanilla Matlab/Python.
25 / 50
Sundials solvers

Python wrappers exist (currently only works for Python ¡= 3.11)

From Matlab2023b, support for CVODE/IDAS solvers available in
Matlab, including sensitivity analysis
Nice blogpost from Mathworks:
https://blogs.mathworks.com/matlab/2024/04/17/
faster-ordinary-differential-equations-odes-solvers-and-se
The above post reports significant performance gains compared to
previous standard ODE solvers in Matlab.

26 / 50
Simulation software (a slightly Julia-centric view)

from https://github.com/SciML/DifferentialEquations.jl
27 / 50
State-space method

y ′ = f (y , z)
0 = g (y , z), gz invertible

Implicit function theorem then gives that there exists (locally) a function
G (y ) such that
z = G (y )
Substitution into the first equations gives the ODE

y ′ = f (y , G (y ))

which can be solved using your method of choice, no new theory. You can
even use an explicit solver if you like.
Lose the structure of the problem which might lead to unnecessary
numerical difficulties.
28 / 50
Outline

Simulation of high index DAE:s, key problems

Simulation of index 1 DAE:s

State-space method
ϵ-embedding
BDF

Index reduction
Index reduction by differentiation
Drift stabilization

29 / 50
ϵ-embedding
Write down, for example, a Runge-Kutta for the ODE
y ′ = f (y , z), ϵz ′ = g (y , z), gz invertible
You then get s
X
Yni = yn + h aij f (Ynj , Znj )
j=1
Xs
ϵZni = ϵzn + h aij g (Ynj , Znj )
j=1
s
X
yn+1 = yn + h bi f (Yni , Zni )
i=1
Xs
ϵzn+1 = ϵzn + h bi g (Yni , Zni )
i=1

30 / 50
ϵ-embedding
Write down, for example, a Runge-Kutta for the ODE

y ′ = f (y , z), ϵz ′ = g (y , z), gz invertible

You then get

(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
ϵ(Zni − zn )i=1,...s = hA(g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
ϵzn+1 = ϵzn + hb T (g (Ynj , Znj ))j=1,...,s

30 / 50
ϵ-embedding
Write down, for example, a Runge-Kutta for the ODE

y ′ = f (y , z), ϵz ′ = g (y , z), gz invertible

You then get

(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
ϵ(Zni − zn )i=1,...s = hA(g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
ϵzn+1 = ϵzn + hb T (g (Ynj , Znj ))j=1,...,s

Assume an implicit method which gives that the A matrix is invertible

h(g (Ynj , Znj ))j=1,...,s = ϵA−1 (Zni − zn )i=1,...s

30 / 50
ϵ-embedding, cont’.
Now, let ϵ → 0 and RK becomes
(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
ϵ(Zni − zn )i=1,...s = hA(g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
ϵzn+1 = ϵzn + hb T (g (Ynj , Znj ))j=1,...,s

31 / 50
ϵ-embedding, cont’.
Now, let ϵ → 0 and RK becomes
(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
ϵ(Zni − zn )i=1,...s = hA(g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
ϵzn+1 = ϵzn + ϵb T A−1 (Zni − zn )i=1,...s

31 / 50
ϵ-embedding, cont’.
Now, let ϵ → 0 and RK becomes
(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
0 = g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
zn+1 = zn + b T A−1 (Zni − zn )i=1,...s

31 / 50
ϵ-embedding, cont’.
Now, let ϵ → 0 and RK becomes
(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
0 = g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
0 = g (yn+1 , zn+1 )

31 / 50
ϵ-embedding, cont’.
Now, let ϵ → 0 and RK becomes
(Yni − yn )i=1,...s = hA(f (Ynj , Znj ))j=1,...,s
0 = g (Ynj , Znj ))j=1,...,s
yn+1 = yn + hb T (f (Ynj , Znj ))j=1,...,s
0 = g (yn+1 , zn+1 )

for a stiffly accurate solver it holds that zn+1 = Zns and the last
rewrite is not necessary.
The solution is identical to the state-space method.
Methods pretty similar but has some pros and cons respectively
State-space form does not require an implicit solver
ϵ-embedding technique possible to generalize to systems not in
semi-explicit form (see Hairer-Wanner)
My ′ = f (t, y )
31 / 50
Order of solver: RK for semi-explicit index 1

If the method is stiffly accurate, the solution is equivalent to

y ′ = f (y , G (y )) and a p order method gives

yn − y (tn ) = O(hp ), zn − z(tn ) = O(hp )

under a Lipschitz assumption on G .

If the solver is not stiffly accurate, you may lose order on the
z-component.
really the same phenomenon, order reduction, as for stiff ODE:s
Example on page 269 in Ascher-Petzold.

32 / 50
Order reduction

Stiffly accurate is sufficient for semi-explicit index 1

Does not apply for higher index
Stiff decay DIRK (Diagonally Implicit Runge Kutta) methods gets
serious order reduction for semi-explicit DAE:s of index 2

33 / 50
Outline

Simulation of high index DAE:s, key problems

Simulation of index 1 DAE:s

State-space method
ϵ-embedding
BDF

Index reduction
Index reduction by differentiation
Drift stabilization

34 / 50
BDF

A s step BDF can be directly applied to the general problem

F (t, y ′ , y ) = 0

without modification with respect to the ODE case.

Popular method. ”BDF is so beautiful that it is hard to imagine
something else could be better”, Petzold, 1988.
IDAS/DASSL(DDASRT/DASPK/DASKR), . . .
Last time I checked OpenModelica used ddasrt (Dubbel precision,
dassl with root solver), a predecessor to DASKR

35 / 50
DASSL

DASSL (and successors) is perhaps the most used DAE solver

Designed to solve DAE:s with index 0 and 1 in the general form

F (t, y ′ , y ) = 0

BDF of order 1 to 5. No order reduction

Variable step length by an extension of fix step length BDF
Will spend time on lecture to describe the basics

36 / 50
Outline

Simulation of high index DAE:s, key problems

Simulation of index 1 DAE:s

State-space method
ϵ-embedding
BDF

Index reduction
Index reduction by differentiation
Drift stabilization

37 / 50
Index reduction

There are many methods to reduce index.

Index reduction through differentiation
Change of variables; think pendulum in polar coordinates. But which
coordinate change? Differential-geometry.
The basic state-space form from control theory

ẋ = f (x, u)
y = h(x, u)

dummy-derivatives, will come back to this next time where automatic

methods suitable for large scale models (Modelica) is discussed.
...

38 / 50
Index reduction by differentiation
Solving high index problems is difficult. For a DAE with index k

F (t, ẏ , y ) = 0

we can derive an ODE by differentiating the equations k times

F (y , ẏ ) = 0
d
F (y , ẏ ) = 0
dt
..
.
dk
F (y , ẏ ) = 0
dt k
This DAE has exactly the same solution set as the original DAE. One
major problem: it is overdetermined!
Find the underlying ODE and simulate that one?
39 / 50
Index reduction by differentiation

Consider the index 1 DAE

ẋ1 = f (x1 , x2 )
0 = g (x1 , x2 ), gx2 (x1 , x2 ) invertible

It is direct to differentiate the second equation and derive the ODE

ẋ1 = f (x1 , x2 )
ẋ2 = −gx2 (x1 , x2 )−1 gx1 (x1 , x2 )f (x1 , x2 )

What will happen with the solution?

40 / 50
Index reduction through differentiation

We have done this before, an index 3 example is

e1 : ẋ = w , e3 : ẏ = z, e5 : 0 = x 2 + y 2 − l 2
e2 : mẇ = −xλ, e4 : mż = −y λ − mg

Differentiate the algebraic equation 2 times and we have an index 1 DAE.

Solution sets to the two
Initial conditions a difficult problem
Underlying ODE (UODE) and the original DAE
Invariants, which are maintained?
Requires projections or other more or less advanced techniques to
fulfill the original algebraic constraints.

41 / 50
Drift in pendulum

0.4 0.4

0.2 0.2

0.0 0.0

−0.2 −0.2
y

y
−0.4 −0.4

−0.6 −0.6

−0.8 −0.8

−1.0 −1.0
−0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6
x x

2.5 sec. 250 sec.

42 / 50
Drift in pendulum
Drift in the pendulum example
Drift error

0.20

0.15
Length error

0.10

0.05

0.00

0 50 100 150 200 250

t [s]

43 / 50
Drift

Theorem
If we apply a method of order p we will (in the example from the last slide)

∥x 2 + y 2 − 1∥ ≤ hp (Atn + Btn2 )

What can you do about this drift?

Baumgarte stabilization
Projection based methods
Use another index reduction technique

44 / 50
Baumgarte stabilization

The first (1976) method to stabilize drift. The principle is simple, instead
of using the second derivative of the algebraic constraint

g̈ = 0

in the solver you use

g̈ + αġ + βg = 0
where α and β are chosen such that the zeros of the polynomial

s 2 + αs + β

lies in the left half plane.

Simple to generalize. Kan be tricky to choose parameters α and β with
respect to stiffness and other numerical properties.

45 / 50
Baumgarte, example

Consider an index 2 DAE with corresponding differentiated index 1 DAE

e1 : ẋ = f (x, y ) e1 : ẋ = f (x, y )
e2 : 0 = g (x) ė2 : 0 = ġ (x) = g ′ (x)f (x, y )

Simulating {e1 , ė2 } will have problems ensuring g (x) = 0. Instead,

simulate the index 1 DAE

ẋ = f (x, y )
0 = g + αġ (x) = g (x) + αg ′ (x)f (x, y )

and then then

g + αg ′ = 0 ⇒ g ∼ e −αt
Thus, the constant α > 0 stabilizes g .

46 / 50
Stabilization through projection

Show the basic principle on a semi-explicit DAE with index 2. Not easy to
generalize for higher index, see Hairer-Wanner for further discussions.

e1 : y ′ = f (y , z)
e2 : 0 = g (y )

Differentiate once
0 = gy (y )f (y , z)
By solving an index 1 DAE (e1 , e2′ ) we will not necessarily fulfill g (yn ) = 0
at each step, even if we start in a consistent starting point.

47 / 50
Stabilization through projection, cont.

Principle
1 Start in a point yn−1 , zn−1 .
2 Take a step to ỹn , z̃n with any method.
3 Project! One projection that has been suggested is defined by

min ∥ỹn − yn ∥, g (yn ) = 0

yn

This is a non-linear optimization’s problem with constraints.

There are many other ways to project to the surfaceM.

48 / 50
ODE:s with invariants

y ′ = f (y ), φ(y ) = 0

Invariants from conservation laws, index reduction

Difference compared to DAE, over determined
φ(y ) = 0 is called a first integral if φ(y )f (y ) ≡ 0 in the neighborhood
of the solution.
Linear first integrals is fulfilled for most methods of integration
Quadratic first integrals is fulfilled by, e.g., symplectic Runge-Kutta
More complex invariants are normally not fulfilled.

49 / 50
 Lecture 2 – Simulation of differential-algebraic equations
Simulation of DAE models and index reduction

Erik Frisk
erik.frisk@liu.se

Department of Electrical Engineering

Linköping University

May 21, 2024

50 / 50