Document [original]

TECHNISCHE UNIVERSIT¨

AT BERLIN

Regularization and Rothe Discretization of

Semi-Explicit Operator DAEs

Robert Altmann Jan Heiland

Preprint 2016/02

Preprint-Reihe des Instituts f¨ur Mathematik

Technische Universit¨at Berlin

http://www.math.tu-berlin.de/preprints

Preprint 2016/02 February 2016

REGULARIZATION AND ROTHE DISCRETIZATION OF

SEMI-EXPLICIT OPERATOR DAES

R. ALTMANN∗, J. HEILAND†

Abstract. A general framework for the regularization of constrained PDEs, also called

operator differential-algebraic equations (operator DAEs), is presented. The given pro-

cedure works for semi-explicit and semi-linear operator DAEs of first order including the

Navier-Stokes and other flow equations. The proposed reformulation is consistent, i.e.,

the solution of the PDE remains untouched. Its main advantage is that it regularizes the

operator DAE in the sense that a semi-discretization in space leads to a DAE of lower

index. Furthermore, a stability analysis is presented for the linear case which shows that

the regularization provides benefits also for the application of the Rothe method. For

this, the influence of perturbations is analyzed for the different formulations. The results

are verified by means of a numerical example with an adaptive space discretization.

Key words. PDAE, operator DAE, regularization, index reduction, Rothe method, method of lines, perturbation

analysis

AMS subject classifications. 65J08,65M12,65L80

1. Introduction

Constrained PDEs arise naturally in the modelling of physical, chemical, and many other

real-world phenomena. They occur whenever different PDE models are coupled, e.g., via

mutual variables at the interfaces, since the coupling is typically modelled via algebraic

constraints. Such models are widely used in flexible multibody dynamics, e.g., the pan-

tograph and catenary benchmark problem [AS00] or the flexible slider crank mechanism

[Sim96, Sim06]. Also flow equations such as the Navier-Stokes equations [Tem77, Wei97]

can be seen as constrained PDEs due to the coupling of momentum equation to the

divergence-free constraint. Further applications can be found in circuit simulation [Tis04],

electromagnetics, and chemical engineering [CM99].

We consider these equation systems of ordinary or partial differential equations (ODEs,

PDEs) and algebraic equations in line with other constrained PDEs – often referred to as

PDAEs – as differential-algebraic equations (DAEs) in function spaces, so-called abstract

or operator DAEs.

Despite the large range of applications and the advantages from the modeling perspec-

tive, the mathematical analysis of operator DAEs is full of open research questions. There

is still no common classification like the index concepts for DAEs [LMT13, Ch. 12]. The

generalization of the tractability index as proposed, e.g., in [Tis04] does not apply for the

commonly used formulation by means of Gelfand triples. The very general concept of the

perturbation index, as it was defined in [RA05] for linear PDAEs, applies under strong reg-

ularity conditions but is still ambiguous in the choice of the norm in which one measures

the perturbation and their derivatives. Also the differentiation index was generalized to

PDAEs [MB00] but has difficulties with the agreement of the PDAE index with the index

Date: February 26, 2016.

The work of R. Altmann was supported by the ERC Advanced Grant ”Modeling, Simulation and

Control of Multi-Physics Systems” MODSIMCONMP and the Berlin Mathematical School BMS.

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 3

of the semi-discretized DAE. Yet another idea is to classify the index of a PDAE directly

by the index that may be determined after a spatial discretization. This, however, leads

to the similar unclear problem, what a good discretization of a PDAE is.

Within this paper, we analyse constrained systems of first order and semi-explicit struc-

ture, particularly we consider systems of the form

˙u(t) + Ku(t) + B∗λ(t) = F(t),(1.1a)

Bu(t) = G(t).(1.1b)

Therein, λdenotes the Lagrange multiplier which enforces the linear constraint Bu=G.

In view of numerical simulations, the incorporation of the constraints via a Lagrange

multiplier and a suitable reformulation seem promising and follows the paradigm in the

treatment of DAEs, that it is preferable to collect all available information in the form of

constraints instead of eliminating them. For the Navier-Stokes equations this means to

maintain the pressure as part of the system.

In this paper, we introduce a regularization or index reduction method for the PDAE

(1.1) without introducing an index as such. We rather refer to the well-defined index

of the semi-discrete system after a spatial discretization by mixed finite elements. In

other words, we propose a reformulation of the given PDAE system such that a semi-

discretization leads to a DAE of lower index, as it was done for second-order systems

appearing in elastodynamics presented in [Alt13]. The advantage is that a transformation

on operator level can be the base for numerically advantageous discretization schemes. The

commonly taken approach of first discretizing and then transforming the equations comes

with the latent risk that the algebraic manipulation are not valid in infinite dimensions

[Hei14]. This may cause instabilities or inconsistencies as the discretization becomes more

accurate.

The main contribution of this paper is then the analysis of the Rothe discretization

through the application of the implicit Euler scheme to the operator DAE (1.1). As for

the finite-dimensional case, we expect a different behavior of the variables uand λand it

will turn out, that we need stronger regularity assumptions to prove the convergence of

the Lagrange multiplier. Among others, we consider the influence of perturbations and

quantify them in a general convergence result for the Rothe method. We show that the

proposed reformulation improves the robustness against such perturbations, as we confirm

numerically for a simulation setup with adaptive, and thus changing, meshes.

The paper is organized as follows. In Section 2, we provide the theoretical framework

and recall the notion of Gelfand triples and Nemytskii maps which are basic tools in

functional analysis for the formulation of operator differential equations. These tools are

then used for the formulation and regularization of the operator DAEs in Section 3. Here

we provide a general framework for linear time-dependent constraints and analyse the

effect of perturbations. The advantages of the obtained formulation and the justification

of calling this procedure an index reduction on operator level is topic of Section 4.

In Section 5 we consider the discretization in time which corresponds to the Rothe

method for time-dependent PDEs. We prove the convergence of the Euler scheme and dis-

cuss the resulting advantages in terms of perturbations. Finally, we illustrate the obtained

theoretical results in numerical simulations of the Navier-Stokes equations in Section 6 and

conclude the paper in Section 7.

4 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

2. Preliminaries

This section is devoted to the introduction of spaces and operators which are needed

for the analysis in Sections 3 and 5 below. Throughout this paper, we use the notion of

Sobolev spaces as in [AF03] and Bochner spaces as in [Rou05, Ch. 1.5].

2.1. Sobolev-Bochner Spaces. To keep the setting as general as possible, we consider

a real, separable, and reflexive Banach space Vand a real separable Hilbert space H

with inner product (·,·). We assume that the spaces V,H, and V∗form a Gelfand triple

(also called evolution triple) [Zei90, Ch. 23.4]. This means that Vis densely, continuously

embedded in H, written as V,→ H, and that Hand its dual space H∗are identified via

the Riesz isomorphism. Such a triple implies the inclusion H∗,→ V∗in the sense that for

h∈ H ∼

=H∗and v∈ V we have

hh, viV∗,V= (h, v).

The space for the Lagrange multiplier is denoted by Qand is assumed to be a real, sep-

arable, and reflexive Banach space. The constraint operator Bthen maps from Vto Q∗.

Together with its dual operator B∗, we obtain the following diagram:

V,→ H =H∗,→ V∗

Q∗Q

BB∗

Example 2.1. A typical example for a Gelfand triple V,→ H ,→ V∗is given by the

Sobolev spaces V:= H1

0(Ω), H:= L2(Ω), and V∗=H−1(Ω).

We consider time derivatives in the generalized sense as defined, e.g., in [Zei90, Ch. 23.5].

We require solutions of system (1.1) to satisfy

u∈Lp(0, T;V) with ˙u∈Lq(0, T;V∗),

where 1 < q ≤p < ∞. If qis the conjugated exponent, i.e., 1/p + 1/q = 1, then, by the

well-known embedding theorems for Gelfand triples [Zei90, Th. 23.23], it holds that such

a solution uis continuous as a function u: [0, T]→ H, i.e., u∈C([0, T],H). Thus, an

initial condition u(0) = afor a∈ H is well-defined.

Remark 2.2.The regularization proposed in Section 3 operates with splittings of the state

space Vand is independent of the time regularity of the function uor ˙u. Thus, we can

also consider less regular systems with ˙u∈Lq(0, T;V∗) with q≤1−1/p, as they may

appear in applications. However, we will have to assume the well-posedness of the initial

condition.

Remark 2.3.Assuming that q≤pis no restriction for applications, where typically p≥2

and, thus, the conjugated exponent is smaller than 2. In this case, from the boundedness

of the interval (0, T) and from the continuous embeddings V,→ H ,→ V∗, one can deduce

that v∈Lp(0, T;V) is also a function in Lq(0, T;V∗).

2.2. Operator K.Consider a possibly nonlinear operator K: (0, T)× V → V∗and let

1≤q,p < ∞. The question arises whether this operator induces a (bounded) operator of

the form

K:Lp(0, T;V)→Lq(0, T;V∗),

(Ku)(t) := K(t, u(t)).

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 5

If such an operator exists, then we do not distinguish between these two notions. We state

a well-known result for Nemytskij mappings for the considered setup of abstract functions.

Theorem 2.4 (cf. [Rou05, Thm. 1.43]).If the operator K: (0, T)×V → V∗is such that

(a) K(t, ·): V → V∗is continuous for almost all t∈(0, T),

(b) K(·, v): (0, T )→ V∗is measurable for all v, and

Vfor some γ∈Lq(0, T),

then the mapping defined via

(Kv)(t) := K(t, v(t)),

is continuous as a map K:Lp(0, T;V)→Lq(0, T;V∗), where 1≤p < ∞and 1≤q≤ ∞.

The case that the exponents 1 < p,q < ∞are conjugated, i.e, 1/p + 1/q = 1, is often

assumed for the analysis of nonlinear evolution equations with monotonicity arguments

[Rou05, Ch. 2 and Ch. 8]. However, for nonlinear operators, even if they are uniformly

bounded as a map V → V∗, the conjugacy of the time exponents may not hold a priori

[Emm04, Ch. 8.2].

Example 2.5 (Navier-Stokes operator).Consider the nonlinear operator which arises in

the weak formulation of the Navier-Stokes equations,

K:V → V∗,hKu, wiV∗,V:= Zu·∇·uw dx.

Then, K:V → V∗is bounded independently of t, cf. [Tem77, Lem. II.1.1], but, in the

three-dimensional case, it is only bounded as an operator K:L2(0, T;V)∩L∞(0, T;H)→

L4/3(0, T;V∗), see e.g. [Rou05, Ch. 8.8.4].

Example 2.6 (p-Laplacian).For the p-Laplacian, i.e.,

Ku, vV∗,V:= ZΩ|∇u|p−2∇u·∇vdx,

we take the Sobolev space V=W1,p

0(Ω). This then induces an operator K:Lp0, T;V)→

Lp0(0, T;V∗) with 1/p + 1/p0= 1, see [Ruˇz04, Ch. 3.3.6].

For special operators K, as, e.g., linear operators that are uniformly bounded with

respect to time, we state the following result.

Corollary 2.7. Consider 1≤p < ∞and an operator K: (0, T)× V → V∗which is

measurable for fixed v∈ V and uniformly bounded in the sense that there exists a constant

CKsuch that kK(t)vkV∗≤CKkvkVfor all v∈ V and almost all t∈(0, T). Then,

(Kv)(t) := K(t, v(t)) defines a continuous operator from Lp(0, T;V)to Lp(0, T;V∗).

Proof. The application of Theorem 2.4 with p=qand γ= 0 yields the result. 

Example 2.8 (Linear elasticity).In the case of linear isotropic material laws for a d-

dimensional domain Ω, i.e.,

Ku, vV∗,V:= ZΩ2µε(u) + λtrace ε(u)Id×d:ε(v) dx

with ε(u) denoting the symmetric gradient, with the Lam´e constants µ,λ, cf. [BS08,

Ch. 11], and A:B:= Pi,j AijBij, we use as ansatz space V=H1(Ω). This setting then

implies an operator K:L2(0, T;V)→L2(0, T;V∗).

6 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

3. Regularization of Operator DAEs

In this section, we consider semi-explicit operator equations in a time interval (0, T ).

With a given linear constraint incorporated by the Lagrangian method, we obtain a system

of the form: find u: (0, T )→ V and λ: (0, T )→ Q such that

˙u(t) + K(t)u(t) + B(t)∗λ(t) = F(t) in V∗,(3.1a)

B(t)u(t) = G(t) in Q∗

(3.1b)

for t∈(0, T) a.e. with initial condition

u(0) = a∈ H.(3.1c)

Therein, B∗denotes the dual of the linear constraint operator B. The operator differential

equation (3.1a) with constraint B(t)u(t) = G(t) is a generalization of a semi-explicit DAE

since here, u(t) belongs to the infinite-dimensional Banach space V. Because of this, we

call system (3.1) a semi-explicit operator DAE.

Suitable function spaces for the solution (u, λ) will be discussed in Theorem 3.8 below.

We assume F ∈ Lq(0, T ;V∗) and G ∈ Lp(0, T;Q∗). The equalities (3.1a) and (3.1b) should

be understood pointwise in L1

loc in the corresponding dual product. By the fundamental

theorem of variational calculus [Emm04, Thm. 8.1.3] and the definition of the weak time

derivative, this means that ˙v(t) = F(t) in V∗if

−ZT

0v(t), wV∗,V˙

φ(t) dt=ZT

0F(t), wV∗,Vφ(t) dt

for all w∈ V and φ∈ C∞

0(0, T). Furthermore, we assume operators K:Lp(0, T;V)→

Lq(0, T;V∗), cf. Section 2.2, and B:Lp(0, T ;V)→Lp(0, T ;Q∗) with 1 < p ≤q < ∞.

3.1. Assumptions on B.In this subsection, we summarize the properties of the con-

straint operator Bwhich we require for a reformulation of the operator DAE (3.1). Note

that we do not need additional assumptions of Kat this point.

Assumption 3.1 (Properties of B).The constraint operator B(t): V → Q∗satisfies the

following conditions:

(a) B(t) is linear and uniformly bounded, and B(·)vis measurable for all v∈ V,

(b) VB:= ker B(t) is independent of time t,

bounded operator B−(t): Q∗→ V such that for all q∈ Q∗it holds that

B(t)B−(t)q=q,

(d) the range of the right inverse Vc:= range B−(t) is independent of time t, and

(e) there exist continuous time derivatives ˙

B(t): V → Q and ˙

B−(t): Q∗→ V.

Remark 3.2 (Time-independent constraint).If the constraint operator is independent of

time, i.e., B(t)≡ B, then Assumption 3.1 reduces to the points (a) and (c).

Remark 3.3 (Induced operators).By Corollary 2.7 it follows that B(t) and B−(t) from

Assumption 3.1 induce bounded operators of the form

B:Lp(0, T;V)→Lp(0, T;Q∗) and B−:Lp(0, T;Q∗)→Lp(0, T;Vc).

Note that the choice of the right inverse in Assumption 3.1 is not unique. A special case,

for which the existence of a right-inverse is guaranteed, is when B(t) satisfies an inf-sup

condition of the form

inf

q∈Q

sup

v∈V

hB(t)v, qi

kvkVkqkQ≥β > 0.

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 7

Nevertheless, this does not imply the time-independence of the range of B−(t). In the next

lemma, we summarize several properties of the right inverse B−(t) from Assumption 3.1.

Lemma 3.4 (Properties of B−).Let Bsatisfy Assumption 3.1. Then, the right inverse

B−(t): Q∗→ V is linear and one-to-one. Furthermore, Vc:= range B−(t)is a closed

subspace of Vand the operator B−(t)B(t): V → V, restricted to Vc, equals the identity.

Proof. The linearity of B−(t) follows from the linearity of the operator B(t) [RR04, Ch. 8.1.2].

For the one-to-one relation, consider q1,q2∈ Q∗with B−(t)q1=B−(t)q2. Then, the ap-

plication of B(t) yields q1=B(t)B−(t)q1=B(t)B−(t)q2=q2.

The linearity of B−(t) and the continuity of B−(t) and B(t) imply that Vcis a closed

subspace of V. Finally, for v∈ Vcand fixed t∈(0, T) there exists q∈ Q∗with B−(t)q=v.

Then, Assumption 3.1 implies

v=B−(t)q=B−(t)B(t)B−(t)q=B−(t)B(t)v. 

Remark 3.5.Lemma 3.4 implies that B−(t)B(t): V → V is a projection onto Vc. Further-

more, the induced Nemytskij mapping B−is a right inverse of B.

An important implication of Assumption 3.1 and Lemma 3.4 is the decomposition of

Lp(0, T;V) as given in the following lemma. This decomposition will be the basis for the

index reduction procedure of Section 3.2.

Lemma 3.6 (Decomposition of Lp(0, T;V)).Consider the subspaces VBand Vcof Vfrom

Assumption 3.1. Then, we have the decomposition

Lp(0, T;V) = Lp(0, T;VB)⊕Lp(0, T;Vc).

Proof. For given v∈Lp(0, T;V), we define r:= Bv∈Lp(0, T ;Q∗), cf. Remark 3.3. Then,

a decomposition of v∈Lp(0, T;V) is given by

v=v0+vc:= v−B−r+B−r.(3.2)

Obviously, vc=B−r∈Lp(0, T;Vc) and v0∈Lp(0, T;VB) follows from Assumption 3.1 by

Bv0=Bv− BB−Bv= 0. We show that the decomposition in (3.2) is unique. For this,

consider v0,w0∈Lp(0, T;VB) and vc,wc∈Lp(0, T;Vc) with v=v0+vc=w0+wc. The

application of Byields Bvc=Bwc. Furthermore, there exist rv,rw∈Lp(0, T;Q∗) such

that vc=B−rvand wc=B−rw. By Assumption 3.1 we then obtain

rv−rw=BB−rv−BB−rw=Bvc−Bwc= 0.

Thus, it holds that vc=B−rv=B−rw=wcand finally also v0=w0.

Within the last lemma of this subsection we consider time derivatives of functions in

Lp(0, T;V).

Lemma 3.7. Let Wbe a closed subspace of Vsuch that there exists a projection P:V → V

that maps Vonto Wand consider v∈Lp(0, T;W). Then, the existence of a time derivative

˙v∈Lp(0, T;V)implies ˙v∈Lp(0, T;W).

Proof. Assume v∈Lp(0, T ;W) with ˙v∈Lp(0, T;V). By assumption, it holds that

(id −P)v(t) = 0 for almost all t∈(0, T) with id denoting the identity. Since the time

derivative of vexists (at least in a generalized sense [Zei90, Ch. 23.5]), we can write

(id −P)˙v(t) = 0, which implies for t∈(0, T ) a.e.,

˙v(t) = P˙v(t)∈ W.

8 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

3.2. Reformulation. This subsection is devoted to the reformulation of the operator

DAE (3.1). In Section 3.3 we discuss the resulting positive effects in terms of the sensitivity

to perturbations. In Section 4 we then show that the presented reformulation is in fact an

index reduction on operator level.

We adapt the technique of minimal extension [KM04], which is an index reduction

procedure especially suitable for semi-explicit DAEs. For this, we first add to system (3.1)

the time derivative of the constraint,

B˙u+˙

Bu=˙

Clearly, this requires the right-hand side Gto be differentiable in the generalized sense,

i.e., G ∈ W1,p(0, T;Q∗). Note that this assumption is already needed for the existence of

a solution of (3.1). This fact comes from the theory of DAEs, see for example [KM06,

Th.2.29] which shows that even for the finite dimensional case with constant coefficients

higher derivatives of the right-hand side are necessary. At this point, also ˙u∈Lp(0, T;V)

seems to be a necessary condition. However, as the next paragraph shows, this requirement

applies only to a part of ˙u.

Second, we use the decomposition from Lemma 3.6 to split uinto u1∈Lp(0, T;VB) and

u2∈Lp(0, T;Vc). Therewith, the two constraints reduce to

Bu2=Gand B˙u2+˙

Bu2=˙

Thus, it is sufficient that the derivative of u2is an element of V. For uin general, we only

need that ˙u∈Lq(0, T;V∗). The assumed regularity of Gimplies with Assumption 3.1,

Lemma 3.4, and equation (3.1b) that u2∈W1,p(0, T;Vc).

Having added one equation, in a third step, we introduce a new variable v2:= ˙u2∈

Lp(0, T;Vc). Recall that Vcis a subspace of Vfor which there exists a projection, cf.

Lemma 3.4. Thus, we can apply Lemma 3.7 at this point. The addition of a new variable

compensates the redundancy of the two constraints. Note that in the reformulated system

the variable u2is not differentiated anymore such that we only need an initial condition

for u1. The initial condition for u2in the original formulation corresponds to a consistency

condition, which typically appears for DAEs [KM06, Ch. 1]. In the sequel, we omit the

time-dependency of the operators Kand B. The overall system then reads: for data

F ∈ Lq(0, T;V∗) and G ∈ W1,p(0, T;Q∗) find functions u1∈Lp(0, T;VB) with ˙u1∈

Lq(0, T;V∗), u2,v2∈Lp(0, T;Vc), and λ∈Lp0(0, T;Q) such that

˙u1(t) + v2(t) + Ku1(t) + u2(t)+B∗λ(t) = F(t) in V∗,(3.3a)

Bu2(t) = G(t) in Q∗,(3.3b)

Bv2(t) + ˙

Bu2(t) = ˙

G(t) in Q∗

(3.3c)

holds for t∈(0, T) a.e. with the initial condition

u1(0) = a0:= a−B−G(0) ∈ H.(3.3d)

The initial condition is well-posed for time differentiable G, since W1,p(0, T ;Q∗) is con-

tinuously embedded in the space of continuous functions with values in Q∗, namely

C([0, T],Q∗) [Rou05, Lem. 7.1]. Further, note that Lp0(0, T;Q) is the right space for

the multiplier λ, since for a separable Banach space Q∗the dual space of Lp0(0, T;Q) can

be identified with Lp(0, T;Q∗), cf. [Rou05, Prop. 1.38].

In the following theorem, we discuss the connection of the original system (3.1) and the

regularized formulation (3.3).

Theorem 3.8 (Equivalence of the reformulation).Consider exponents 1< q ≤p < ∞,

and p0with 1/p0+ 1/p = 1. Assume that F ∈ Lq(0, T;V∗),G ∈ W1,p(0, T ;Q∗), and

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 9

a∈ H as well as the operator Bsatisfying Assumption 3.1. Then, the operator DAE (3.1)

has a solution (u, λ)with u∈Lp(0, T;V),˙u∈Lq(0, T;V∗), and λ∈Lp0(0, T ;Q)if and

only if system (3.3) has a solution (u1, u2, v2, λ)with u1∈Lp(0, T, VB),˙u1∈Lq(0, T;V∗),

u2,v2∈Lp(0, T;Vc), and λ∈Lp0(0, T;Q). Furthermore, it holds that u=u1+u2and

˙u2=v2.

Proof. Let (u, λ) be a solution of (3.1). We define

u1:= u−B−Bu∈Lp(0, T;VB) and u2:= B−Bu∈Lp(0, T;Vc).

With equation (3.1b), we obtain u2=B−Gand thus, by the regularity of Gand Assump-

tion 3.1, ˙u2∈Lp(0, T;Vc). With v2:= ˙u2the quadruple (u1, u2, v2, λ) satisfies equations

(3.3a-c). The initial condition (3.3d) is satisfied because of

u1(0) = u(0) −u2(0) = a−B−G(0).

For the reverse direction consider a solution of (3.3), namely (u1, u2, v2, λ). Then,

u:= u1+u2∈Lp(0, T;V) and because of the regularity of Gand equation (3.3b), by

Remark 2.3, it holds that ˙u= ˙u1+ ˙u2∈Lq(0, T;V∗). We show that ˙u2=v2. Equation

(3.3c) and the time derivative of equation (3.3b) yield

Bv2+˙

Bu2=˙

G=d

dtBu2=B˙u2+˙

Bu2.

Note that ˙u2∈Lp(0, T ;Vc), as shown in the first part of the proof. The invertibility of B

on Vc, see Lemma 3.4, then gives ˙u2=v2. Thus, the pair (u, λ) satisfies equations (3.1a)

and (3.1b). For the initial condition (3.1c), we obtain

u(0) = u1(0) + u2(0) = a−B−G(0) + B−G(0) = a. 

From the solution representation given in Theorem 3.8 we deduce that not every initial

condition a∈ H admits a solution to (3.1).

Corollary 3.9. Let the assumptions of Theorem 3.8 hold. For the existence of a solution

to (3.1) it is necessary that the initial data a∈ H can be decomposed as a=a0+B−G(0),

where B−G(0) ∈ Vcand a0is in the closure of VBin H.

We discuss two examples for which we obtain different kinds of consistency conditions.

Example 3.10. If the operator Bequals the divergence operator and V= [H1

0(Ω)]d, then

VBdenotes the space of divergence-free functions in V. In this case, the closure of VBwith

respect to H= [L2(Ω)]dis a proper subspace of H, cf. [Tem77, Ch. 1, Thm. 1.4],

H=v∈ H | ∇·v= 0, v ·ν∂Ω= 06=H.

Therein, ν∂Ωdenotes the normal outer vector along the boundary. Note that the closure

is even a subspace of H(div,Ω) = {v∈ H | ∇ · v∈L2(Ω)}. Thus, the initial value a0

cannot be chosen arbitrarily in H.

Example 3.11. If Bequals the trace operator, i.e., B:V:= H1(Ω) →H1/2(∂Ω), then we

have VB=H1

0(Ω). Since the closure of H1

0(Ω) in H:= L2(Ω) equals Hitself, the initial

data only has to satisfy a0∈ H. This means, however, that the initial data a=a0+B−G(0)

can also be chosen arbitrarily in Hsuch that there is no consistency condition.

10 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

3.3. Influence of Perturbations. As mentioned in the introduction, we do not define an

index for operator DAEs. Nevertheless, the influence of perturbations provides information

about the stability of the system, similar to the perturbation index for DAEs.

In this subsection, we show the positive effect of the presented regularization in terms

of perturbations. For this, we restrict the analysis to the case p=q= 2 with a linear,

symmetric, on VBelliptic, and bounded operator K, i.e, for u∈ VBand v,w∈ V we

assume

hKu, ui ≥ k1kuk2and hKv, wi ≤ k2kvkkwk.

Note that we use k·k := k·kVand later |·| := k·kHto simplify the notation. We consider

the to (3.1) corresponding perturbed problem

¯u+K¯u+B∗¯

λ=F+δin V∗,(3.4a)

B¯u=G+θin Q∗.(3.4b)

Here, (¯u, ¯

λ) denotes the solution if we include perturbations δ: [0, T]→ V∗and θ: [0, T]→

Q∗. For the regularized equations, the perturbed problem has the form

ˆu1+ ˆv2+K(ˆu1+ ˆu2) + B∗ˆ

λ=F+δin V∗,(3.5a)

Bˆu2=G+θin Q∗,(3.5b)

Bˆv2+˙

Bˆu2=˙

G+ξin Q∗.(3.5c)

For this system, we consider perturbations of the form δ∈L2(0, T;V∗) and θ,ξ∈

L2(0, T;Q∗) and the solution is denoted by (ˆu1,ˆu2,ˆv2,ˆ

λ). The initial condition is given

by ˆu1(0) = u1(0) −e1,0, i.e., e1,0contains the initial error. Because of Theorem 3.8, it

is sufficient to consider the regularized system (3.5). The result for the original operator

DAE then follows if we replace ξby ˙

θ, cf. Remark 3.13 below.

By Cemb we denote the continuity constant of the embedding V,→ H. Furthermore, we

introduce the errors

e1:= ˆu1−u1, e2:= ˆu2−u2, ev:= ˆv2−v2, eλ:= ˆ

λ−λ.

Theorem 3.12. Consider the perturbed problem (3.5) with a linear, symmetric, on VB

elliptic, and bounded operator K. Furthermore, let Bsatisfy Assumption 3.1 and the

perturbations δ∈L2(0, T ;V∗)and θ,ξ∈L2(0, T;Q∗). Then, the error of the differential

variable satisfies with a positive constant c∈Rthat

ke1k2

C([0,T];H)+k1ke1k2

L2(0,T;V)≤ |e1,0|2+chkδk2

L2(0,T;V∗)+kθk2

L2(0,T;Q∗)+kξk2

L2(0,T;Q∗)i.

Proof. Bounds of the errors e2and evcan be easily found by the continuity of ˙

Band

the right-inverse of B. For the error in the differential part u1, we test the difference of

equations (3.5a) and (3.3a) by e1∈ VB. Thus, the term with the Lagrange multiplier

vanishes and the properties of the operator Kcan be exploited. The details can be found

in [Alt15, Sect. 6.1]. 

Remark 3.13.In order to transfer these results to the perturbation analysis of the original

formulation (3.1) we have to insert ξ=˙

θ. Thus, the error also depends on the derivative

of the perturbation θwhich leads to possible instabilities known from high-index DAEs.

If the perturbation is not smooth enough, it may even make the model useless.

In the given setting of the evolution equations, it is not possible to gain similar estimates

for eλ. Estimates of the error in the Lagrange multiplier are only possible if we consider

the primitive of eλor assume more regular perturbations δ∈L2(0, T;H∗) and e1,0∈ VB.

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 11

4. Spatial Discretization

For the simulation of time-dependent PDEs, we need discretizations in time and space.

Because of the special role of the time variable in DAEs we do not consider space-time

schemes. Instead, we consider the approach of discretizing in space and time separately.

In this section, we consider the systems which result from a spatial discretization of the

operator DAE. Thus, we follow the method of lines [Hol07, Ch. 3.4]. The Rothe method

[Rot30], in which one discretizes in time first, is then discussed in Section 5.

As mentioned above, a spatial discretization of an operator DAE leads to a classical

DAE, for which the differentiation index is well-defined [BCP96, HW96, KM06]. Within

this paper, we simply write index, meaning the differentiation index. The index quanti-

fies the necessary number of differentiation steps, in order to obtain an ODE and thus,

describes to which degree the solution depends on derivatives of the involved quantities.

The dependence on derivatives may lead to instabilities within the numerical simulation.

For a precise definition, we refer to [BCP96, Def. 2.2.2]. Recall that we do not use any

index definition of PDAEs.

We show that the DAE corresponding to the original system (3.1) is of index 2 whereas

the DAEs resulting from the reformulated systems are of index 1. For this, only standard

assumptions on the used finite element schemes have to be considered.

4.1. Finite Element Discretization. For the spatial discretization, we consider finite

dimensional approximations of the spaces VB,Vc, and Q. We denote the approximation

spaces by VB,h,Vc

h, and Qh, respectively. Furthermore, we define Vh=VB,h ⊕Vc

has finite

dimensional approximation of V.

Thinking of finite elements on a regular mesh Tof the domain Ω, cf. [Cia78, Bra07],

we consider basis functions {ϕi}1,...,n1of VB,h,{ϕi}n1+1,...,n of Vc

h, and {ψi}1,...,m of Qh

with m=n−n1. Hence, we assume that dim Vc

h= dim Qh. The finite dimensional

approximations of u1, u2, v2, and λare then represented by the coefficient vectors q1, q2, r2,

and µ, respectively. By q∈Rnwe denote the vector q= [qT

1, qT

2]T. Based on this

discretization scheme, we define the positive definite mass matrix M∈Rn,n by Mi,j :=

(ϕi, ϕj). The discrete version of the constraint operator Bis defined by

B(t)∈Rm,n, Bj,i(t) := B(t)ϕi, ψj.(4.1)

Note that, according to Assumption 3.1, it is natural to assume that Bis continuously

differentiable with respect to time and that Bhas full rank.

Remark 4.1 (Nonconforming discretization).In order that Bis well-defined, the operator

Bhas to be defined for the given basis functions. Since nonconforming finite elements are

not excluded [BS08, Ch. 10], the application of Bmay be generalized to a piecewise (with

respect to the triangulation T) application of the operator.

Remark 4.2 (Inf-sup stability).For the unique solvability of the semi-discrete systems

resulting from the finite element discretization it is sufficient that the constraint matrix

Bis of full rank. For a stable approximation of the Lagrange multiplier λin terms of

the discretization parameter h, one may assume an inf-sup condition, i.e., there exists a

constant βdisc >0, independent of hand time t, such that

inf

qh∈Qh

sup

vh∈Vh

hB(t)vh, qhi

kvhk kqhkQ≥βdisc,

cf., e.g., [Bra07, Ch. III.4].

12 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

Finally, we define the discrete version of the operator Kby

K:Rn→Rn,K(q), ek:= K(Xn

i=1qiϕi), ϕk.

Here, ekdenotes the k-th unit vector in Rn. As before, Kcan be written as a n×nmatrix

if Kis linear.

4.2. Discretization of (3.1).For the computation of the index of the resulting DAE

system, we assume that the mass matrix Mis positive definite and that Bis of full rank.

First, we consider the index of the DAE which results from a spatial discretization of the

original system (3.1). With the introduced notation, the DAE has the form

M˙q+K(q) + BT(t)µ=f,(4.2a)

B(t)q=g.(4.2b)

From [HW96, Ch. VII.1] we infer that the system (4.2) is of index 2, since BM−1BTis

invertible. The index-2 structure can also made visible through a differentiation of the

constraint (4.2b). This then leads to the (analytically) equivalent DAE

M BT(t)

B(t) 0 ˙q

µ=f−K(q)

˙g−˙

B(t)q.

Again the assumptions on Mand Bimply that the matrix on the left-hand side is invert-

ible. Thus, a single differentiation leads to an ODE for qand an algebraic equation for

the Lagrange multiplier µ.

4.3. Discretization of (3.3).In the case of a conforming discretization, i.e., VB,h ⊂ VB,

h⊂ Vc, and Qh⊂ Q, the matrix B(t) has the special structure B(t) = [0 B2(t)].

Therein, the matrix B2(t) is square and non-singular. In this case, the semi-discrete

version of (3.3) has the form

M˙q1

r2+K(q1, q2) + 0

2(t)µ=f,(4.3a)

B2(t)q2=g,(4.3b)

B2(t)r2= ˙g−˙

B2(t)q2.(4.3c)

At this point we only mention that this system forms a DAE of index 1. We give the proof

in Lemma 4.3 below for the more general case.

In many cases, one depends on a nonconforming spatial discretization [BS08, Ch. 10],

i.e., the discrete ansatz spaces are not subspaces of the original search spaces. One sim-

ple example is the Crouzeix-Raviart element [CR73], a lowest order piecewise linear but

discontinuous discretization scheme. Since we do not assume VB,h ⊂ VBfor general mixed

finite element discretizations, i.e., ker B6⊂ ker B, cf. [GR86, Ch. 3], we lose the special

structure of B(t).

In general, we have B(t) = [B1(t)B2(t)] and simply assume that the block B2is

non-singular. This is no restriction, since one may always permute the columns of B

(corresponds to a reordering of basis functions in VB,h and Vc

h) such that the B2block is

regular. Then, the semi-discretized system reads

M˙q1

r2+K(q1, q2) + BT(t)µ=f,(4.4a)

B2(t)q2=g−B1(t)q1,(4.4b)

B2(t)r2= ˙g−B1(t) ˙q1−˙

B1(t)q1−˙

B2(t)q2.(4.4c)

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 13

Lemma 4.3 (Index-1 DAE).For a positive definite mass matrix Mand a continuously

differentiable constraint matrix Bwith a non-singular block B2, the DAEs (4.3) and (4.4)

are of index 1.

Proof. Similar to the proof of [KM06, Th. 6.12], we show that (4.4) is of index 1. The

property then follows for system (4.3) as well because it is a special case.

Since the matrix B2(t) is of full rank, equations (4.4b) and (4.4c) yield direct expressions

of q2and r2in terms of q1and ˙q1. Furthermore, a multiplication of (4.4a) from the left

by BM−1provides a formula for µin terms of q1. Here we use the assumptions on M

and Bwhich imply that the matrix BM−1BTis invertible. Finally, inserting all these

expressions into equation (4.4a), we obtain an ODE in q1. Thus, we can solve system (4.4)

without any further differentiation steps. 

4.4. Commutativity. The connection between the presented regularization on operator

level in Section 3 and the index reduction for DAEs is illustrated in Figure 4.1, cf. [Alt13]

for second-order systems. The scheme shows that regularization and spatial discretiza-

tion commute if corresponding discretization schemes are used and the index reduction is

performed by minimal extension. Note that this is beneficial for adaptive simulations as

modifications of the finite element schemes such as a refinement of the underlying mesh

do not call for another index reduction step afterwards. This fact is also of importance

for the Rothe discretization as shown in the numerical example in Section 6.

operator DAE of

’index-2 type’ (3.1)

operator DAE of

’index-1 type’ (3.3)

DAE of index 2 (4.2) DAE of index 1 (4.3)

regularization

of operator DAEs

index reduction by

minimal extension [KM04]

spatial

discretization

spatial

discretization

Figure 4.1. Illustrative scheme of the commutativity of regularization and discretization.

Since the index of the DAE (4.3) is, compared to the DAE (4.2), reduced by one, we may

call the proposed regularization procedure from Section 3 an index reduction on operator

level.

5. Temporal Discretization

For the time integration of the operator DAE (3.3) we restrict ourselves to the implicit

Euler method. We prove the convergence of the scheme and highlight the needed ad-

justments in contrast to operator ODEs (time-dependent PDEs without constraint) for

which the convergence is well-known. Again we benefit from the regularization introduced

in Section 3. Furthermore, we consider perturbations of the right-hand sides in order to

analyse the convergence of the Rothe method applied to operator DAEs.

14 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

5.1. Time-discrete Systems. As in Section 3.3, we restrict the analysis to the linear

case with p=q= 2, i.e., we consider a linear operator Kwhich is symmetric, continuous,

and positive on VB. Furthermore, we assume Bto be time-independent with right inverse

B−:Q∗→ Vcand continuity constant CB−, cf. Assumption 3.1. Recall that (·,·) := (·,·)H

denotes the inner product in H.

Remark 5.1.The given assumptions on the operator Kalready imply that VBis a Hilbert

space. In many applications (·,·) + hK·,·i defines an inner product in V. In this case,

we may assume that the splitting V=VB⊕ Vcis orthogonal with respect to this inner

product. This then allows to prove stronger convergence results.

For the temporal discretization we consider an equidistant partition 0 = t0< t1<··· <

tn=Tof the interval [0, T] with time step size τ. The semi-discrete approximations of

u1,u2,v2, and λat time tj=jτ are then denoted by uj

1,uj

2,vj

2, and λj, respectively.

For the application of the implicit Euler scheme to system (3.3) we replace the derivative

˙u1by the discrete derivative Duj

1:= (uj

1−uj−1

1)/τ. This then leads to stationary PDEs

which have to be solved in every time step. The differential equation (3.3a) turns into

(Duj

1, v) + hKuj

1, vi+hλj,Bvi=hFj, vi−(vj

2, v)−hKuj

2, vi(5.1)

for j= 1, . . . , n, whereas the constraints (3.3b) and (3.3c) result in

hBuj

2, qi=hGj, qi,hBvj

2, qi=h˙

Gj, qi.(5.2)

Therein, we assume that uj−1

1∈ H is given and search for uj

1∈ VB,uj

2, vj

2∈ Vc, and

λj∈ Q. The test functions are given by v∈ V and q∈ Q.

Since Gis continuous, we set Gj=G(tj). Note, however, that Fjand ˙

Gjcannot equal

the function evaluations of Fand ˙

Gat time tj, since this is not well-defined. Instead, we

use integral means over a time interval, i.e.,

Fj:= 1

τZtj

tj−1F(t) dt∈ V∗,˙

Gj:= 1

τZtj

tj−1

G(t) dt∈ Q∗.

With these approximations we define Fτ,Gτ, and ˙

Gτas the piecewise constants

Fτ(t) := Fj,Gτ(t) := Gj,˙

Gτ(t) := ˙

Gj,

for t∈]tj, tj+1] and continuous extension in t= 0. This then implies

Fτ→ F in L2(0, T;V∗),Gτ→ G,˙

Gτ→˙

Gin L2(0, T;Q∗).(5.3)

Since the operator Bis invertible on Vc, we argue from (5.2) that uj

2=B−Gjand

2=B−˙

Gj. Inserting this into (5.1) and testing only with functions v∈ VB, we obtain

(Duj

1, v) + hKuj

1, vi=hFj, vi−(B−˙

Gj, v)−hKB−Gj, vi.(5.4)

This equation will be used for the stability estimates in the following subsection. It

remains to discuss the solvability of system (5.1) for uj

1and λj. With the bilinear form

c:V ×V → R, given by c(u, v) := τ−1(u, v) + hKu, vi,and the functional F∈ V∗,

hF, vi:= hFj, vi−(B−˙

Gj, v)−hKB−Gj, vi+1

τ(uj−1

1, v),

equation (5.4) can be written as c(uj

1, v) = F(v) for all v∈ VB. The unique solvability of

(5.4) for uj

1then follows by the Lax-Milgram lemma [Eva98, Sect. 6.2.1]. For the unique

solvability of λjwe consider equation (5.1) tested by functions v∈ Vc,

hB∗λj, vi=hFj, vi−(B−˙

Gj, v)−hKB−Gj, vi−(Duj

1, v)−hKuj

1, vi.(5.5)

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 15

Equation (5.4) implies that the right-hand side of (5.5) vanishes for all functions in VB.

Thus, the functional is an element of the polar set Vo

Bon which the operator B∗is invertible

[Bra07, Ch. III, Lem. 4.2].

5.2. Stability Estimates. The most important ingredient for the convergence analysis

of the Rothe method are stability or a priori estimates which we provide in this subsection.

Since equation (5.4) is essentially an operator ODE for uj

1, the shown bounds follow the

lines of the stability results in [Emm01, Ch. 4], see also [Tem77, Ch. III.4]. Amongst

others, we take advantage of the equality

2(Duj, uj) = D|uj|2+τ|Duj|2.(5.6)

This identity follows by a simple calculation, cf. [Emm01, Lem. 3.2.2]. For the differential

variable uj

1we obtain the following result.

Lemma 5.2 (Stability).Assume F ∈ L2(0, T;V∗

B)and G ∈ H1(0, T;Q∗). Then, the

approximations uk

1∈ VBgiven by the Euler scheme (5.4) with u0

1∈ H satisfy for all

1≤k≤nthe estimate

|uk

1|2+τ2

j=1 |Duj

1|2+τk1

j=1 kuj

1k2≤M2

(5.7)

with the constant M2:= |u0

1|2+ 3kFk2

L2(0,T;V∗

B)+C2

B−(C4

emb +k2

2)kGk2

H1(0,T;Q∗)/k1.

Proof. Using as test function v=uj

1∈ VB,j≥1, in the Euler scheme (5.4), we obtain

(Duj

1, uj

1) + hKuj

1, uj

1i=hFj, uj

1i−(B−˙

Gj, uj

1)−hKB−Gj, uj

1i.(5.8)

Summation over j= 1, . . . ,k, together with property (5.6), the Cauchy-Schwarz inequal-

ity, and the continuous embedding V,→ H yields

|uk

1|2−|u0

1|2+τ2

j=1|Duj

1|2+ 2τk1

j=1 kuj

1k2

(5.6)

≤2τ

j=1

(Duj

1, uj

1)+2τ

j=1hKuj

1, uj

(5.8)

≤2τ

j=1 kFjkV∗

B+Cemb|B−˙

Gj|+k2kB−Gjkkuj

1k.

The application of Young’s inequality [Eva98, App. B] and the boundedness of B−shows

|uk

1|2−|u0

1|2+τ2

j=1|Duj

1|2+τk1

j=1 kuj

1k2

≤3τ

j=1 kFjk2

V∗

B+C2

B−C4

embk˙

Gjk2

Q∗+C2

B−k2

2kGjk2

Q∗.

Finally, the assertion follows by properties of the Bochner integral which imply that

τPk

j=1 kFjk2

V∗

B≤ kFk2

L2(0,T;V∗

B).

16 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

With the same assumptions as in Lemma 5.2, we may also prove that there exists a

positive constant c∈Rsuch that

j=1 kDuj

1k2

V∗

B≤cM2.(5.9)

This result follows from equation (5.4) which yields for j≥1,

kDuj

1kV∗

B:= sup

v∈VB,kvk=1 hFj, vi−(B−˙

Gj, v)−hKB−Gj, vi−hKuj

1, vi

≤ kFjkV∗

B+CB−C2

embk˙

GjkQ∗+k2CB−kGjkQ∗+k2kuj

1k.

An application of Young’s inequality, the summation over j, and the estimate (5.7) finally

imply (5.9).

Remark 5.3.Assume that (·,·)+hK·,·i defines an inner product in Vwith respect to which

the decomposition V=VB⊕Vcis orthogonal. If we assume more regularity of the given

data in the form of F ∈ L2(0, T;H∗) and u0

1∈ VB, then we obtain the estimate

j=1 |Duj

1|2≤M2

reg.

Here, the constant Mreg includes the terms ku0

1kas well as kFkL2(0,T ;H∗). For this, equation

(5.4) has to be tested by Duj

1∈ VB. Thus, we obtain the estimate (5.9) in a stronger norm

which is crucial in view of the Lagrange multiplier.

In order to prove the convergence of the Euler scheme, we need to define global approx-

imations. Given uj

1,j= 0, . . . , n, we define the piecewise constant and piecewise linear

functions U1,τ and ˆ

U1,τ on the interval [0, T] by

U1,τ (t) := (u0

1if t= 0,

1if t∈]tj−1, tj],ˆ

U1,τ (t) := (u0

1if t= 0,

1+ (t−tj)Duj

1if t∈]tj−1, tj].(5.10)

We show that the sequences U1,τ and ˆ

U1,τ are uniformly bounded. The boundedness of

both sequences in L∞(0, T;H) and L2(0, T;VB) follows directly from (5.7) if we assume

1∈ VB. Additionally, we obtain by the stability estimate (5.9) that



˙

U1,τ 



L2(0,T;V∗

B)=

j=1 Ztj

tj−1kDuj

1k2

V∗

Bdt=τ

j=1 kDuj

1k2

V∗

B≤cM2.(5.11)

For the approximations of u2and v2we define similarly the piecewise constant functions

U2,τ (t) := uj

2if t∈]tj−1, tj], V2,τ (t) := vj

2if t∈]tj−1, tj](5.12)

with a continuous extension in t= 0. Note that this definition implies that U2,τ =

B−Gτand V2,τ =B−˙

Gτ. Finally, we define as approximation of the Lagrange multiplier

Λτ: ]0, T]→ Q by

Λτ(t) := λjif t∈]tj−1, tj].(5.13)

With this, we can rewrite equation (5.5) in the form

hB∗Λτ, vi=hFτ, vi−(B−˙

Gτ, v)−hKB−Gτ, vi−(˙

U1,τ , v)−hKU1,τ , vi.(5.14)

Here, we consider test functions v∈ V and t∈(0, T) a.e..

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 17

Since the Lagrange multiplier corresponds to the algebraic variable in the finite-dimen-

sional case, we expect less regularity than for the other variables. Indeed, we are not able

to bound Λτin L2(0, T;Q) uniformly within the given weak setting. Note that with the

additional regularity assumptions from Remark 5.3, we would be able to show the desired

boundedness. Instead, we consider the primitive of Λτwhich then leads to a weaker notion

of solvability, cf. [EM13]. We define ˜

Λτ∈AC([0, T]; Q) by

Λτ(t) := Zt

Λτ(s) ds.(5.15)

With an integration of (5.14) over [0, t], we obtain an equation for the primitive of the

Lagrange multiplier,

hB∗˜

Λτ, vi=h˜

Fτ, vi−(B−˜

Gτ, v)−hKB−˜

Gτ, vi−(ˆ

U1,τ , v)−hK˜

U1,τ , vi+ (u0

1, v).(5.16)

Therein, ˜

Fτ,˜

Gτ,˜

Gτ, and ˜

U1,τ denote the primitives of Fτ,Gτ,˙

Gτ, and U1,τ , respectively.

For ˜

Λτwe are able to prove the boundedness independently of the step size τ.

Lemma 5.4 (Boundedness of ˜

Λτ).Assume F ∈ L2(0, T ;V∗),G ∈ H1(0, T;Q∗), and

1∈ VB. Then, the sequence ˜

Λτis bounded in C([0, T ]; Q).

Proof. We make use of the inf-sup condition of the operator Band, by equation (5.16),

we obtain the estimate

βinfk˜

ΛτkC([0,T];Q)≤max

t∈[0,T]

sup

v∈V

hB∗˜

Λτ(t), vi

kvk

(5.16)

≤max

t∈[0,T]hk˜

Fτ(t)kV∗+Cemb|B−˜

Gτ(t)|+k2kB−˜

Gτ(t)k

+Cemb|ˆ

U1,τ (t)|+k2k˜

U1,τ (t)k+Cemb|u0

1|i.

By the properties of the Bochner integral and (5.7), this is bounded uniformly in terms of

T, the initial data, and the right-hand sides. For details we refer to [Alt15, Sect. 10.3]. 

5.3. Passing to the Limit. Since every bounded sequence has a weakly convergent sub-

sequence, the results from the previous subsection imply the existence of weak limits U1,

U2,V2, and ˜

Λ. More precisely, we obtain by (5.3) that

U2,τ =B−Gτ→U2:= B−G, V2,τ =B−˙

Gτ→V2:= B−˙

Gin L2(0, T;Vc).

The embedding H1(0, T;Q∗),→C([0, T]; Q∗) implies additionally that U2satisfies the

consistency condition U2(0) = B−G(0). By the boundedness of U1,τ and ˆ

U1,τ , we obtain

weak limits in L2(0, T;VB). Furthermore, the estimate (5.7) implies that



U1,τ −ˆ

U1,τ 



L2(0,T;H)≤τ

j=1 uj

1−uj−1

1

2≤τM2→0.

Thus, the two limits coincide in L2(0, T;H) and the continuous embedding V,→ H implies

that the same is true for the limit in L2(0, T;VB). We denote the joined limit by U1, i.e.,

U1,τ ,ˆ

U1,τ * U1in L2(0, T;VB).

Finally, Lemma 5.4 implies the existence of a weak limit ˜

Λ for which

Λτ*˜

Λ in Lp(0, T;Q)

18 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

for all 1 <p<∞. It remains to show that the obtained limits solve the operator DAE

(3.3). For this, we assume u0

1=a0∈ VB. Obviously, the limits U2and V2solve equations

(3.3b) and (3.3c). For U1we obtain the following result.

Theorem 5.5. Assume F ∈ L2(0, T;V∗

B),G ∈ H1(0, T ;Q∗), and u0

1=a0∈ VB. Then,

the weak limit U1∈L2(0, T;VB)of the sequence U1,τ solves equation (3.3a) in V∗

B, i.e.,

for test functions in VB. Furthermore, U1has a generalized derivative which satisfies

U1∈L2(0, T;V∗

B).

Proof. As in [Emm01, Ch. 4], we consider equation (5.4) with test functions v∈ VBin the

form

dt(ˆ

U1,τ , v) + hKU1,τ , vi=hFτ, vi−(B−˙

Gτ, v)−hKB−Gτ, vi.

With Φ ∈C∞

0(0, T), we may rewrite this in integral form, i.e.,

0−ˆ

U1,τ , v˙

Φ(t) + KU1,τ , vΦ(t) dt=ZT

0Fτ−B−˙

Gτ−KB−Gτ, vΦ(t) dt.

We advance to the limit τ→0 and obtain for the right-hand side

0Fτ−B−˙

Gτ−KB−Gτ, vΦ(t) dt−→ ZT

0F −V2−KU2, vΦ(t) dt.

Furthermore, the weak convergence of U1,τ and ˆ

U1,τ in L2(0, T;VB) implies

0−ˆ

U1,τ , v˙

Φ + KU1,τ , vΦ dt−→ ZT

0−U1, v˙

Φ + KU1, vΦ dt.

As a result, the obtained limit U1∈L2(0, T;VB) satisfies for all v∈ VB,

dt(U1, v)+(U2, v) + hK(U1+U2), vi=hF, vi.(5.17)

It remains to show that U1has a generalized derivative. From the definition of ˆ

U1,τ in

(5.10) we know that its time derivative equals Duj

1for t∈]tj−1, tj[ and is bounded in

L2(0, T;V∗

B) due to (5.11). Thus, there exists a subsequence which weakly converges to

a limit V1∈L2(0, T;V∗

B). For every Φ ∈C∞

0(0, T) and v∈ VBthis limit satisfies the

equality

0U1(t), v˙

Φ(t) dt= lim

τ→0ZT

0ˆ

U1,τ (t), v˙

Φ(t) dt

= lim

τ→0−ZT

0˙

U1,τ (t), vΦ(t) dt=−ZT

0V1(t), vΦ(t) dt.

This shows that ˙

U1=V1∈L2(0, T;V∗

B) in the generalized sense. Finally, we have to check

whether U1satisfies the stated initial condition. Since ˆ

U1,τ * U1as well as d

dtˆ

U1,τ *˙

U1=

V1, for Φ ∈C1([0, T]) with Φ(T) = 0 and arbitrary v∈ VB, we derive by the integration

by parts formula that

0 = lim

τ→0ZT

0˙

U1,τ −˙

U1, vΦ dt=−a0−U1(0), vΦ(0).

Since VBis dense in VB

Hby definition, this implies U1(0) = a0.

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 19

Since we were not able to show the uniform boundedness of Λτ, we cannot prove its

convergence to the weak solution of the operator DAE. However, we can show that the

limit ˜

Λ solves the operator DAE in a weaker sense. To be more precise, we state this result

in the following theorem.

Theorem 5.6. Assume F ∈ L2(0, T;V∗),G ∈ H1(0, T ;Q∗), and u0

1=a0∈ VB. Then, for

any sequence of step sizes with τ→0the sequence ˜

Λτconverges weakly to ˜

Λin L2(0, T;Q)

and (U1, U2, V2,˜

Λ) solves system (3.3) in the weak distributional sense, meaning that for

all v∈ V and Φ∈C∞

0(0, T)it holds that

0−U1, v˙

Φ + V2, vΦ + K(U1+U2), vΦ−B∗˜

Λ, v˙

Φdt=ZT

0F, vΦdt.

Remark 5.7.If we assume more regularity of the data as in Remark 5.3, i.e., F ∈

L2(0, T;H∗) and the orthogonality of the decomposition V=VB⊕ Vcwith respect to

the inner product (·,·) + hK·,·i, then the (weak) limits U1,U2,V2, and Λ solve the regu-

larized operator DAE (3.3). In addition, we have ˙

U1∈L2(0, T;H).

5.4. Perturbations. In the previous subsection, we have proven the convergence of the

Euler scheme if we assume that the PDEs were solved exactly in every time step. In order

to prove the convergence of the Rothe method, errors due to the spatial discretization have

to be included as well. For this, we consider the time-discrete systems with additional

perturbations of the right-hand sides. This may then be interpreted as the error of a

spatial discretization, cf. [Alt15, Sect. 10.4].

We consider system (3.3), discretized by the implicit Euler scheme, with additional per-

turbations. Note that we still consider the linear case, cf. Section 5.1 for the assumptions

on Kand B. The differences of the exact and perturbed solution (ˆuj

1,ˆuj

2,ˆvj

2,ˆ

λj), namely,

1:= ˆuj

1−uj

1∈ VB, ej

2:= ˆuj

2−uj

2∈ Vc, ej

v:= ˆvj

2−vj

2∈ Vc, ej

λ:= ˆ

λj−λj∈ Q

then satisfy the equations

Dej

1+ej

v+Kej

1+ej

2+B∗ej

λ=δjin V∗,(5.18a)

Bej

2=θjin Q∗,(5.18b)

Bej

v=ξjin Q∗.(5.18c)

As in the previous section, the reachable results depend on the assumed smoothness of

the given data. Since we want to include estimates for the Lagrange multiplier, which

underlines the positive effects of the regularization from Section 3, we consider the more

regular case. For this, we consider perturbations δj∈ H∗and θj,ξj∈ Q∗. Note that

for an estimate of the errors ej

1,ej

2, and ej

vit would be sufficient to assume δj∈ V∗.

Furthermore, we assume the spaces VBand Vcto be orthogonal with respect to the inner

product defined by (·,·) + hK·,·i.

By equations (5.18b) and (5.18c) we directly obtain the estimates

kej

2k ≤ CB−kθjkQ∗,kej

vk ≤ CB−kξjkQ∗.(5.19)

For estimates on ej

1, we may test equation (5.18a) by ej

1, similarly as in Lemma 5.2.

Because of the additional smoothness δ∈ H∗we may also test equation (5.18a) by Dej

and obtain, due to the assumed orthogonality of VBand Vc,

|Dej

1|2+hKej

1, Dej

1i=hδj, Dej

1i−(ej

v, Dej

1)+(ej

2, Dej

1).

20 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

The equality 2hKej

1, Dej

1i=DhKej

1, ej

1i+τhKDej

1, Dej

1ithen yields together with Young’s

inequality,

|Dej

1|2+DhKej

1, ej

1i+τk1kDej

1k2≤3kδjk2

H∗+ 3C2

embkej

vk2+ 3C2

embkej

2k2.

A summation for j= 1, . . . , k and a multiplication by τfinally leads to

k1kek

1k2+τ

j=1|Dej

1|2≤k2ke0

1k2+ 3τ

j=1 kδjk2

H∗+C2

embkej

2k2+C2

embkej

vk2.(5.20)

Furthermore, equation (5.18a) and the inf-sup condition of Byield

βinfkej

λkQ≤sup

v∈Vc

hB∗ej

λ, vi

kvk≤ kδjkV∗+k2kej

1k+k2kej

2k+C2

embkej

vk+Cemb|Dej

1|.

Putting the estimates (5.19) and (5.20) together, we obtain with a generic constant, which

we express by ., that

τβ2

inf

j=1 kej

λk2

(5.19)

.τ

j=1 kδjk2

V∗+kθjk2

Q∗+kξjk2

Q∗+τ

j=1 kej

1k2+|Dej

1|2

(5.20)

.ke0

1k2+τ

j=1 kδjk2

H∗+kθjk2

Q∗+kξjk2

Q∗.(5.21)

To summarize the result for the Lagrange multiplier, we assume that all the perturbations

are of the same order of magnitude, i.e., δj≈δ∈ H∗,θj≈θ∈ Q∗, and ξj≈ξ∈ Q∗. Then,

the piecewise constant function Eλ: [0, T]→ Q, defined by Eλ(t) = ej

λfor t∈]tj−1, tj],

satisfies

βinf kEλkL2(0,T;Q).ke0

1k+√TkδkH∗+kθkQ∗+kξkQ∗.

This estimate raises hope that numerical simulations can be performed in a reasonable

fashion. Note that this is only true for the regularized operator DAE. For an analogous

estimate for the original formulation, equation (5.18c) has to be replaced by BDej

2=Dθj.

Thus, the perturbation ξjhas to be replaced by the discrete derivative of θjwhich then

leads to an additional term that scales with τ−1in the error estimates.

6. Examples

This section is devoted to illustrate the benefits of the regularization investigated in

Section 3 by means of numerical approximations of the Navier-Stokes equations. First, we

revisit the numerical results presented in [AH15] in the context of a Rothe discretization.

Second, we give a concrete example of the perturbation results from Section 5.4 by illus-

trating how a perturbation induced by a mesh adaption gets numerically differentiated in

the index-2 but not in the index-1 formulation.

6.1. Navier-Stokes Equations. We consider the standard formulation of the Navier-

Stokes equations [Tem77] for an incompressible flow in a domain Ω ⊂Rd,

(6.1) ˙u+ (u·∇)u−ν∆u+∇p=f, ∇·u= 0.

We interpret the pressure pas a Lagrange multiplier that couples the incompressibility

constraint ∇·u= 0 to the state equations. Then, equation (6.1) takes the form of system

(3.1) with the spaces chosen as V= [H1

0(Ω)]d,H= [L2(Ω)]d, and Q=L2(Ω)/R.

The operator B:V → Q∗is defined as the divergence and K:V → V∗is the operator

representing convection and diffusion, cf. Example 2.5. Without further assumptions, the

nonlinearity Konly extends to K:L2(0, T;V)→L1(0, T;V∗), cf. [Tem77, Lem. III.3.1].

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 21

This causes the main difficulties in the existence theory for the Navier-Stokes equations.

However, this lower regularity does not affect the splitting as proposed in Section 3. In

particular, Assumption 3.1 is fulfilled. The weak divergence operator is linear and bounded

and there exists a continuous right inverse as shown, e.g., in [Tar06, Lem. I.4.1]. The

splitting V=VB⊕ Vcis then given by the space of divergence-free functions and its

(orthogonal) complement. Note that our approach is different from [EM13] where the

splitting of Vis used to eliminate the constraints rather then to augment the system.

For a fixed spatial discretization and for similar approaches to the nonlinear parts, the

discrete equations obtained via Rothe’s method coincide with the equations stemming

from the method of lines. Thus, the numerical study conducted in [AH15] also serves as

an example for the advantages of the index-1 formulation as the base for Rothe’s method.

6.2. Adaptive Changes of the Mesh. The benefits of the index-1 formulation become

particularly apparent for space discretizations that change with time, when, e.g., the mesh

is adapted to the current state of the system. Note that the opportunity to adapt the

mesh between time steps is the major advantage of Rothe’s method over the method of

lines.

Let the superscripts +, c, and −denote the next, current, and previous value of the

variables, respectively. We use the same superscripts for the discrete operators to denote

possibly different spatial discretizations. We consider the algebraic systems that are ob-

tained from (6.1) after the time-discretization and, subsequently, the discretization of the

space on the currently considered mesh. With the same notation as used in [AH15, Ch.

3.3] for the method of lines, for the original system (3.1), the update to (q+, p+) from the

current iterate (qc, pc) via a time step of length τis obtained via

(6.2) 1

τM+−B+T

B+0q+

p+=1

τM+qc+f+−K+(q+)

g+.

To advance by one time-step and the regularized index-1-type formulation (3.3), we pro-

pose the solution of







τM+

11 M+

12 −B+

τM+

21 M+

22 −B+

τB+

1B+

20 0

B+

10 0 B+













˜q+







=





τM+

11qc

1+f+

1−K+

1(q+

1, q+

τM+

21qc

1+f+

2−K+

2(q+

1, q+

τB+

1qc

1+ ˙g+







.(6.3)

The different stability properties become evident, if one examines the inherent equation

for the pressure update p+, derived via premultiplying the upper part of the equations by

B+(M+)−1. In the index-2 case (6.2), this leads to

(6.4) −B+(M+)−1B+Tp+=B+qc−B+q+

τ+B+(M+)−1f+−K(q+).

The index-1 formulation yields for the pressure

(6.5) −B+(M+)−1B+Tp+=1

τB+qc

1−q+

−τ˜q+

2+B+(M+)−1f+−K(q+

1, q+

2).

Comparing equations (6.4) and (6.5) to the time-continuous formula for the pressure

−B+(M+)−1B+Tp+=−˙g++B+(M+)−1f+−K(q+),

we find that the consistency errors are given as

ind2 := −1

τ(B+qc−g+)−˙g+and e+

ind1 := −1

τB+qc

1−q+

−τ˜q+

2−˙g+

22 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

for the index-2 and index-1 scheme, respectively. For different meshes B+qc6=Bcqc=gc,

and, unless the changes in the discretization and in Bare smooth, the consistency error

ind2 will not approach zero as τ→0. More likely, for switches in the mesh, which may

appear in every time step, this term leads to an error in p+that scales with τ−1, cf.

Section 5.4. In the index-1 formulation, the error term e+

ind1 is not present at all since the

equation 1

τB+qc

1−q+

−τ˜q+

2−˙g+= 0 is an explicit part of the numerical scheme.

Remark 6.1.The error e+is due to the changing meshes. In the case of inexact system

solves it will add to the errors that were investigated in [AH15].

For a numerical example, we consider the following numerical approach to the cylinder

wake at Reynolds number Re = 60 as described in [AH15]. The spatial component is

discretized by means of Crouzeix-Raviart elements on a coarse grid with 7404 velocity

and 2413 pressure nodes and on a fine grid with 15110 velocity and 4963 pressure nodes.

The time evolution on the interval [0,2] is discretized using the Euler scheme – implicit in

the linear parts and explicit in the nonlinearity – on a grid of 2048 and 4096 equidistant

points and starting with the steady state Stokes solution. As reference solution we use

the trajectory on the finer spatial mesh and a time discretization by means of the implicit

trapezoidal rule.

To illustrate the error in the pressure induced by changes in the mesh and its scaling

with the inverse of the time step length, we start the simulation on the fine grid and switch

to the coarse grid at t= 0.67. At t= 1.33 we switch back to the fine grid. The code for

this numerical example is available from the author’s public git repository [Hei16].

The error e+

ind2 and its scaling are well visible in the index-2 formulation while the

pressure error at the mesh switches in the index-1 formulation is less prominent and

obviously independent of the time-step size, cf. Figure 6.1(a). Therein we have plotted

the pointwise in time approximation errors kpref(t)−pNtskL2(Ω), where pref is the reference

solution, pNts is the numerical approximation, and Ω is the computational domain. Note,

that in the index-2 case, the error e+

ind2 affects the pressure update only instantaneously.

This is due to the implicit decoupling of pressure and velocity that makes the velocity

approximation independent of the pressure so that an error in the pressure will not spread

in time. Accordingly, the large amplitude in the pressure error at the switching points is

not seen in the velocity approximations, cf. Figure 6.1(b), where the pointwise temporal

error in the velocity, that is defined in the same way as the pressure error, is plotted. For

inexact solves, however, the implicit splitting of qand pin (6.2) is not exact such that

a single occurrence of e+will spread to the velocity computation and, thus, linger in the

velocity approximation forever, cf. the numerical results in [AH15].

Note that the error due to the time approximation is only visible in the beginning. After

the switch to the coarse mesh the spatial error dominates.

7. Conclusion

Within this paper, we have introduced a reformulation for a special class of semi-explicit

operator DAEs with linear constraints such that a spatial discretization by finite elements

leads to a DAE of index 1, rather than index 2. Thus, the procedure can be seen as an

index reduction for operator DAEs.

Furthermore, we have proven the convergence of the Rothe discretization on the base

of the implicit Euler scheme. We have quantified the advantages of the reformulation

over the original schemes in terms of stability estimates concerning the robustness against

perturbations in the right-hand sides. Particularly, we have shown that derivatives of

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 23

0 0.67 1.33 2

10−4

10−3

10−2

10−1

100

time t

(a): pointwise error in the pressure p

Nts = 2048 (index-1)

Nts = 4096 (index-1)

Nts = 2048 (index-2)

Nts = 4096 (index-2)

0 0.67 1.33 2

10−5

10−4

10−3

10−2

time t

(b): pointwise error in the velocity q

Nts = 2048 (index-1)

Nts = 4096 (index-1)

Nts = 2048 (index-2)

Nts = 4096 (index-2)

Figure 6.1. Pointwise in time error for qand pfor the index-1 and index-2

formulations for various number of time steps Nts and with mesh switches

at t= 0.67 and t= 1.33 plotted at every 50-th point of the temporal grid

and at the points where the switches occur.

perturbations, that may occur in the original formulation, are not present in the solutions

of the reformulated equations.

Finally, we have illustrated the advantages of the regularized formulation in a numer-

ical simulation of flow equations where the spatial discretization changes at certain time

points.

References

[AF03] R. A. Adams and J. J. F. Fournier. Sobolev Spaces. Elsevier, Amsterdam, second

edition, 2003.

[AH15] R. Altmann and J. Heiland. Finite element decomposition and minimal exten-

sion for flow equations. ESAIM Math. Model. Numer. Anal., 49(5):1489–1509,

2015.

24 REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES

[Alt13] R. Altmann. Index reduction for operator differential-algebraic equations in

elastodynamics. Z. Angew. Math. Mech. (ZAMM), 93(9):648–664, 2013.

[Alt15] R. Altmann. Regularization and Simulation of Constrained Partial Differential

Equations. PhD thesis, Technische Universit¨at Berlin, 2015.

[AS00] M. Arnold and B. Simeon. Pantograph and catenary dynamics: A benchmark

problem and its numerical solution. Appl. Numer. Math., 34(4):345–362, 2000.

[BCP96] K.E. Brenan, S.L. Campbell, and L. R. Petzold. Numerical solution of initial-

value problems in differential-algebraic equations. Society for Industrial and

Applied Mathematics (SIAM), Philadelphia, PA, 1996.

[Bra07] D. Braess. Finite Elements - Theory, Fast Solvers, and Applications in Solid

Mechanics. Cambridge University Press, New York, third edition, 2007.

[BS08] S. C. Brenner and L. R. Scott. The Mathematical Theory of Finite Element

Methods. Springer-Verlag, New York, third edition, 2008.

[Cia78] P. G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland,

Amsterdam, 1978.

[CM99] S. L. Campbell and W. Marszalek. The index of an infinite-dimensional implicit

system. Math. Comput. Model. Dyn. Syst., 5(1):18–42, 1999.

[CR73] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element

methods for solving the stationary Stokes equations. I. Rev. Franc. Automat.

Inform. Rech Operat, 7(R-3):33–75, 1973.

[EM13] E. Emmrich and V. Mehrmann. Operator differential-algebraic equations arising

in fluid dynamics. Comp. Methods Appl. Math., 13(4):443–470, 2013.

[Emm01] E. Emmrich. Analysis von Zeitdiskretisierungen des inkompressiblen Navier–

Stokes-Problems. Cuvillier Verlag, G¨ottingen, Germany, 2001.

[Emm04] E. Emmrich. Gew¨ohnliche und Operator-Differentialgleichungen. Vieweg, Wies-

baden, Germany, 2004.

[Eva98] L. C. Evans. Partial Differential Equations. American Mathematical Society

(AMS), Providence, second edition, 1998.

[GR86] V. Girault and P.-A. Raviart. Finite Element Methods for Navier–Stokes Equa-

tions. Theory and Algorithms. Springer, Berlin, Germany, 1986.

[Hei14] J. Heiland. Decoupling and optimization of differential-algebraic equations with

application in flow control. PhD thesis, Technische Universit¨at Berlin, 2014.

[Hei16] J. Heiland. TayHoodMinExtForFlowEqns. Public Git Repository, 2016. Solution

of time-dependent 2D nonviscous flow with nonconforming minimal extension,

https://github.com/highlando/TayHoodMinExtForFlowEqns.

[Hol07] M. H. Holmes. Introduction to Numerical Methods in Differential Equations.

Springer-Verlag, New York, 2007.

[HW96] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and

Differential-Algebraic Problems. Springer-Verlag, Berlin, second edition, 1996.

[KM04] P. Kunkel and V. Mehrmann. Index reduction for differential-algebraic equations

by minimal extension. Z. Angew. Math. Mech., 84(9):579–597, 2004.

[KM06] P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and

Numerical Solution. European Mathematical Society Publishing House, Z¨urich,

Switzerland, 2006.

[LMT13] R. Lamour, R. M¨arz, and C. Tischendorf. Differential-Algebraic Equations: A

Projector Based Analysis. Differential-Algebraic Equations Forum. Springer,

Heidelberg, Germany, 2013.

[MB00] W. S. Martinson and P. I. Barton. A differentiation index for partial differential-

algebraic equations. SIAM J. Sci. Comput., 21(6):2295–2315, 2000.

REGULARIZATION AND ROTHE DISCRETIZATION OF SEMI-EXPLICIT OPERATOR DAES 25

[RA05] J. Rang and L. Angermann. Perturbation index of linear partial differential-

algebraic equations. Appl. Numer. Math., 53(2-4):437–456, 2005.

[Rot30] E. Rothe. Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindi-

mensionaler Randwertaufgaben. Math. Ann., 102(1):650–670, 1930.

[Rou05] T. Roub´ıˇcek. Nonlinear Partial Differential Equations with Applications.

Birkh¨auser, Basel, Switzerland, 2005.

[RR04] M. Renardy and R. C. Rogers. An Introduction to Partial Differential Equations.

Springer-Verlag, New York, second edition, 2004.

[Ruˇz04] M. Ruˇziˇcka. Nichtlineare Funktionalanalysis: Eine Einf¨uhrung. Springer-Verlag,

London, 2004.

[Sim96] B. Simeon. Modelling a flexible slider crank mechanism by a mixed system of

DAEs and PDEs. Math. Comp. Model. Dyn., 2:1–18, 1996.

[Sim06] B. Simeon. On Lagrange multipliers in flexible multibody dynamics. Comput.

Methods Appl. Mech. Engrg., 195(50–51):6993–7005, 2006.

[Tar06] L. Tartar. An Introduction to Navier–Stokes Equation and Oceanography.

Springer, New York, NY, 2006.

[Tem77] R. Temam. Navier–Stokes Equations. Theory and Numerical Analysis. North-

Holland, Amsterdam, Netherlands, 1977.

[Tis04] C. Tischendorf. Coupled systems of differential algebraic and partial differential

equations in circuit and device simulation. Modeling and numerical analysis.

Habilitationsschrift, Technische Universit¨at Berlin, 2004.

[Wei97] J. Weickert. Applications of the theory of differential-algebraic equations to par-

tial differential equations of fluid dynamics. PhD thesis, Fakult¨at f¨ur Mathe-

matik, Technische Universit¨at Chemnitz, 1997.

[Zei90] E. Zeidler. Nonlinear functional analysis and its applications. II/A: Linear

monotone operators. Springer, Berlin, Germany, 1990.

∗Institut f¨

ur Mathematik MA4-5, Technische Universit¨

at Berlin, Straße des 17. Juni

136, 10623 Berlin, Germany

†Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstraße 1,

39106 Magdeburg, Germany

E-mail address:[email protected], [email protected]