DAE examples

The following examples of differential algebraic equations are drawn from the book
Peter Kunkel and Volker Mehrmann: Differential Algebraic Equations. Analysis and numerical methods. Zürich: European Mathematical Society Publishing House. (2006).

1 Example 1: A chemical reaction

<a name="test1"></a>

A chemical reaction is controled by external cooling.

c

is the concentration of the reactant,

T

temperature,

R

the reaction rate per unit volume,

T_{0}

the initial temparature,

c_{0}

the initial concentration and

R_{0}

the initial rate.

T_{C}

denotes the external cooling temperature (a given function of time). The equation reads

(\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} \overset{\cdot}{c} \\ \overset{\cdot}{T} \\ \overset{\cdot}{R} \end{matrix}) = (\begin{matrix} k_{1} (c_{0} - c) - R \\ k_{1} (T_{0} - T) + k_{2} R - k_{3} (T - T_{C}) \\ R - k_{3} \exp (- k_{4} / T) c \end{matrix})

k_{i}, i = 1, 2, 3, 4

are given constants which we choose as

k_{1} = 0.9, k_{2} = 1.0, k_{3} = 0.5 and k_{4} = 0.2 .

The standard initial values are

c_{0} = 1, T_{0} = 100 .

For consistency it then follows that

R_{0} = 0.499000999333666 .

R_{0}

is obtained from the third equation using

c_{0}

and

T_{0}

. As

T_{C}

we take

T_{C} = 0

\overset{\cdot}{c}

and

\overset{\cdot}{T}

for

t = t_{0}

are taken from the system, and

\overset{\cdot}{R}

for

t = t_{0}

is obtained by differentiating the third equation w.r.t.

t

and inserting these values. This equation has index one: differentiate the third equation and insert the derivatives from the first two in order to obtain an explicit system. We take the rather long time interval

[0, 5]

. You may change the values

c_{0}, T_{0}, t_{end}

2 Example 2

<a name="test2"></a>

The van der Pol equation

\begin{matrix} \overset{\cdot}{y} & = & z \\ ϵ \overset{\cdot}{z} & = & (1 - y^{2}) z - y \end{matrix}

and its limit case for

ϵ = 0

\begin{matrix} \overset{\cdot}{y} & = & z \\ 0 & = & (1 - y^{2}) z - y \end{matrix}

You may work with smaller and smaller

ϵ

and with

ϵ = 0

. Then this equation has also index 1: differentiate the second equation and insert

z

for

\overset{\cdot}{y}

. We use as standard initial values

y_{0} = - 1 / 10

and

z_{0} = - 100 / 99

. This implies

{\overset{\cdot}{y}}_{0} = - 100 / 99

{\overset{\cdot}{z}}_{0}

is chosen consistently as

0

. For these values and

ϵ = 0

the equation exhibits a singularity very near to

1.805

. We take

1.805

as final time. For

ϵ > 0

there is no such singularity, but cusps build up for smaller and smaller values. You may change here the initial value

y_{0}

and

z_{0}

ϵ > 0

or only

y_{0} \not = \pm 1

ϵ = 0

. In the latter case it follows from the consistency conditions that

\begin{matrix} z_{0} & = & y_{0} / (1 - y_{0}^{2}) \\ {\overset{\cdot}{y}}_{0} & = & z_{0} \\ {\overset{\cdot}{z}}_{0} & = & {\overset{\cdot}{y}}_{0} (1 + 2 y_{0} z_{0}) / (1 - y_{0}^{2}) \end{matrix}

3 Example 3

<a name="test3"></a>

The pendulum, a typical from the constrained motion examples (like occur in multibody motions applications):

\begin{matrix} m \overset{\cdot\cdot}{x} + 2 x λ & = & 0 \\ m \overset{\cdot\cdot}{y} + 2 y λ + mg & = & 0 \\ x^{2} + y^{2} - L^{2} & = & 0 . \end{matrix}

Here

m

is the mass mounted at the tip of the pendulum,

x, y

is the position of the tip and

λ

is the force acting along the pendulum.

g

is the earth acceleration. This equation has index three. It is rewritten in first order form using

\begin{matrix} z_{1} & = & x \\ z_{2} & = & y \\ z_{5} & = & λ \end{matrix}

\begin{matrix} z_{3} - \overset{\cdot}{z_{1}} & = & 0 \\ z_{4} - \overset{\cdot}{z_{2}} & = & 0 \\ m \overset{\cdot}{z_{3}} + 2 z_{1} z_{5} & = & 0 \\ m \overset{\cdot}{z_{4}} + 2 z_{2} z_{5} + mg & = & 0 \\ z_{1}^{2} + z_{2}^{2} - L^{2} & = & 0 . \end{matrix}

In this form it cannot be solved using the BDF-code, since the Jacobian does not allow the computation of

{\overset{\cdot}{z}}_{5}

. We proceed in transforming it into index one:
Repeated differentiation of the fifth equation gives

\begin{matrix} \overset{\cdot}{z_{1}} z_{1} + \overset{\cdot}{z_{2}} z_{2} & = & 0 \\ \Rightarrow z_{3} z_{1} + z_{4} z_{2} & = & 0 \\ \Rightarrow \overset{\cdot}{z_{3}} z_{1} + z_{3} \overset{\cdot}{z_{1}} + \overset{\cdot}{z_{4}} z_{2} + z_{4} \overset{\cdot}{z_{2}} & = & 0 \\ \Rightarrow (- 2 z_{1} z_{5} / m) z_{1} + (- g - 2 z_{2} z_{5} / m) z_{2} + z_{3} \overset{\cdot}{z_{1}} + z_{4} \overset{\cdot}{z_{2}} & = & 0 \\ \Rightarrow L^{2} z_{5} + {mgz}_{2} / 2 - m (z_{3} \overset{\cdot}{z_{1}} + z_{4} \overset{\cdot}{z_{2}}) / 2 & = & 0 \\ \Rightarrow L^{2} \overset{\cdot}{z_{5}} + mg \overset{\cdot}{z_{2}} / 2 - m (\overset{\cdot}{z_{3}} \overset{\cdot}{z_{1}} + z_{3} \overset{\cdot\cdot}{z_{1}} + \overset{\cdot}{z_{4}} \overset{\cdot}{z_{2}} + z_{4} \overset{\cdot\cdot}{z_{2}}) / 2 & = & 0 \\ \Rightarrow L^{2} \overset{\cdot}{z_{5}} + mg \overset{\cdot}{z_{2}} / 2 - m (\overset{\cdot}{z_{3}} z_{3} + z_{4} \overset{\cdot}{z_{4}}) & = & 0 \\ \Rightarrow L^{2} \overset{\cdot}{z_{5}} + mg \overset{\cdot}{z_{2}} / 2 - (z_{3} (- 2 z_{1} z_{5}) + z_{4} (- 2 z_{2} z_{5} - mg)) & = & 0 \\ \Rightarrow L^{2} \overset{\cdot}{z_{5}} + (3 / 2) mg \overset{\cdot}{z_{2}} & = & 0 \end{matrix}

The last equation gives us the required explicit representation of

\overset{\cdot}{z_{5}}

. As final semiexplicit system of index one we now use

\begin{matrix} z_{3} - \overset{\cdot}{z_{1}} & = & 0 \\ m \overset{\cdot}{z_{3}} + 2 z_{1} z_{5} & = & 0 \\ z_{1}^{2} + z_{2}^{2} - L^{2} & = & 0 \\ z_{3} z_{1} + z_{4} z_{2} & = & 0 \\ L^{2} z_{5} + (mg / 2) z_{2} - (m / 2) (z_{3}^{2} + z_{4}^{2}) & = & 0 \end{matrix}

For the computation of consistent initial values you can chose two parameters only:

ϑ_{0}

and

α

. Then

\begin{matrix} x_{0} & = & \cos (ϑ_{0}) = z_{1} \\ y_{0} & = & \sin (ϑ_{0}) = z_{2} \\ {\overset{\cdot}{x}}_{0} & = & \sin (ϑ_{0}) α = {\overset{\cdot}{z}}_{1} = z_{3} \\ {\overset{\cdot}{y}}_{0} & = & - \cos (ϑ_{0}) α = {\overset{\cdot}{z}}_{2} = z_{4} \\ z_{5} & = & (- {mgz}_{2} / 2 + m (z_{3}^{2} + z_{4}^{2}) / 2) / L^{2} \\ {\overset{\cdot}{z}}_{3} & = & - 2 z_{1} z_{5} / m \\ {\overset{\cdot}{z}}_{4} & = & - g - 2 z_{2} z_{5} / m \\ {\overset{\cdot}{z}}_{5} & = & - (3 / 2) {mgz}_{4} / L^{2} \end{matrix}

We use

ϑ_{0} = 5 π / 4

and

α = 0.5

as standard values.

L = 1

and

m = 1

File translated from T_EX by T_TM Unregistered, version 4.03.
On 16 Jun 2016, 18:43.