Document [original]

TECHNISCHE UNIVERSIT¨

AT BERLIN

Convergence of the Rothe Method Applied to

Operator DAEs Arising in Elastodynamics

Robert Altmann

Preprint 2015/20

Preprint-Reihe des Instituts f¨ur Mathematik

Technische Universit¨at Berlin

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

Preprint 2015/20 August 2015

CONVERGENCE OF THE ROTHE METHOD APPLIED TO

OPERATOR DAES ARISING IN ELASTODYNAMICS

R. ALTMANN∗

Abstract. The dynamics of elastic media, constrained by Dirichlet boundary condi-

tions, can be modeled as operator DAE of semi-explicit structure. These models include

flexible multibody systems as well as applications with boundary control. In order to use

adaptive methods in space, we analyse the properties of the Rothe method concerning

stability and convergence for this kind of systems.

For this, we consider a regularization of the operator DAE and prove the weak con-

vergence of the implicit Euler scheme. Furthermore, we consider perturbations in the

semi-discrete systems which correspond to additional errors such as spatial discretization

errors.

Key words. PDAE, operator DAE, regularization, evolution equations, elastodynamics, Rothe method, Euler method

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

1. Introduction

Within this paper, we prove the convergence of the implicit Euler method applied

to a differential-algebraic equation (DAE) in the abstract setting, i.e., to an operator

DAE. More precisely, we consider the dynamics of elastic media which are constrained by

Dirichlet boundary conditions or by a coupling condition. This also includes the simulation

of flexible multibody systems.

The modeling of mechanical systems often leads to constrained systems of ordinary

and partial differential equations (PDEs), see [Sim98, SGS06]. In this paper, we restrict

the analysis to elastodynamics with constraints on the boundary in the form of Dirichlet

boundary conditions. These constraints are incorporated by the Lagrangian method, since

the Dirichlet data may be unknown a priori [Sim06]. Such systems arise typically in the

field of flexible multibody dynamics [Sha97, GC01, Bau10, Sim13]. Therein, a number of

deformable bodies are coupled through joints. Note, however, that the considered approach

also includes a more general coupling of flexible bodies with any other dynamical system

as well as boundary control. This is possible, since the Dirichlet boundary conditions are

explicitly given in the system equations in form of a constraint. A spatial discretization of

such systems then typically leads to DAEs of differentiation index 3 [Sim06]. The concept

of the differentiation index measures, loosely speaking, how far the DAE is apart from

an ordinary differential equation and thus, provides a measure of difficulty. A general

introduction to DAEs can be found in the monographs [CM99, KM06, Ria08, LMT13],

see also the review on the different index concepts in [Meh13].

If we want to analyse a dynamical system before the spatial discretization, the frame-

work of classical DAEs is too restrictive. Since the elastic behaviour is described by PDEs,

we obtain so-called partial differential-algebraic equations (PDAEs) or, formulated in a

weak functional analytic setting, operator DAEs. Following [Zei90, Ch. 23], we use for this

Date: September 3, 2015.

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

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

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 3

formulation appropriate Sobolev-Bochner spaces and different spaces for the solution and

their derivatives.

Although operator DAEs provide a quick and simple modeling procedure, there is no

general theory for the well-posedness or a classification as for DAEs [Tis03, LMT13]. For

the analysis of such systems we need to explore the interaction of DAE theory and operator

theory.

The basis for the presented convergence analysis of the time integration is the regularized

operator formulation introduced in [Alt13]. Therein, the operator DAE describing the

motion of elastic media was reformulated such that a spatial discretization leads to a DAE

of index 1 rather than index 3. The regularization also improves the sensitivity in terms of

perturbations which is of importance for numerical simulations. This then allows to apply

the Rothe method - which is popular for time-dependent PDEs - also for operator DAEs.

This enables adaptive procedures, especially in the space variable, since the underlying

grid may be changed easily from time step to time step.

The paper is organized as follows. In Section 2 we recall the equations of motion

for the dynamics of elastic media and their formulation as operator DAE. Furthermore,

we present the regularized formulation of this system, which corresponds to the index-1

formulation in finite dimensions. The discretization of the time variable is then subject of

Section 3. Therein, we discuss the solvability of the semi-discrete system and provide

stability estimates. To show the convergence of the Rothe method, we define global

approximations and analyze their behaviour in the limit when the step size goes to zero. In

Section 4 we employ these results and consider additional perturbations. This is important

for the convergence of the Rothe method, since in practice additional errors due to a spatial

discretization appear. Finally, we give some concluding remarks in Section 5.

2. Preliminaries

Since we focus on Dirichlet boundary conditions, the trace operator which extends the

mapping

γ:C(¯

Ω) →C(∂Ω), u 7→ γ(u) := u|∂Ω

(2.1)

to Sobolev spaces in of special importance, cf. [AF03]. We introduce the necessary Sobolev

and Sobolev-Bochner spaces of abstract functions. All these spaces are necessary for the

weak formulation of PDEs and their interpretation as operator differential equations. With

this, we formulate the equations of motion as operator DAE in Section 2.3.

2.1. Spaces and Norms. Within this paper, we use the standard notation of Sobolev

spaces, i.e., L2(Ω) denotes the Lebesgue space of square integrable functions and H1(Ω)

the functions with an additional derivative (in the weak sense). Furthermore, H1

ΓD(Ω)

denotes the subspace of H1(Ω) with vanishing trace on ΓD⊆∂Ω, i.e., with homogeneous

boundary values on ΓD, and H1/2(ΓD) the space of traces, cf. [AF03].

For the weak formulation of elastodynamics we need functions in several components.

For this, we define

V:= H1(Ω)d,VB:= H1

ΓD(Ω)d,H:= L2(Ω)d,Q∗:= [H1/2(ΓD)]d.

Note that the Hilbert space Qis only defined via its dual space Q∗. For the inner product

in Hand the norms in Hand Vwe use the abbreviations

(u, v) := (u, v)H,|u|:= kukH,kuk:= kukV.

4 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

Note that the spaces V,H,V∗form a Gelfand triple [Wlo87, Ch. 17.1]. Thus, we have

,−→ H ∼

=H∗d

,−→ V∗

with continuous and dense embeddings. The equivalence of the Hilbert spaces Hand H∗

follows from the Riesz representation theorem [RR04, Th. 6.52]. Within this setting, the

duality pairing of Vand V∗is the continuous extension of the inner product in H, i.e., for

h∈ H and v∈ V, we obtain

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

As solutions of the operator equations below, we consider abstract functions, i.e., func-

tions on a time interval [0, T] which map to a Banach space X. In our case, these Banach

spaces are again the Sobolev spaces introduced above. This leads to Sobolev-Bochner

spaces, see [Rou05, Ch. 7] for an introduction.

In particular, we consider the space L2(0, T;X) which includes abstract functions

u: [0, T]→Xwith

kuk2

L2(0,T;X)=ZT

ku(t)k2

Xdt < ∞.

Similarly, Hk(0, T ;X) denotes the Sobolev-Bochner space of abstract functions which have

time derivatives (in the distributional sense) up to order kin L2(0, T ;X).

The space of continuous functions with values in Xis donoted by C([0, T]; X).

2.2. Elastodynamics. We review the governing equations for the dynamics of elastic

media. Throughout this paper, Ω ⊂Rddenotes a domain with Lipschitz boundary where

ΓD⊆∂Ω denotes the Dirichlet boundary and ΓN=∂Ω\ΓDthe Neumann boundary.

Note that we do not consider the pure Neumann problem, i.e., ΓN=∂Ω, since this would

exclude the considered coupling throughout the boundary.

2.2.1. Equations of Motion. The equations of elastodynamics describe the evolution of a

deformable body under the influence of applied forces based on Cauchy’s theorem [Cia88,

Ch. 2]. We consider the theory of linear elasticity for homogeneous and isotropic materials,

i.e., we assume small deformations only. Note that for large deformations the nonlinear

theory has to be applied in order to obtain reasonable results. The corresponding initial-

boundary value problem in the classical form with prescribed Dirichlet data uDand applied

forces βand τreads

ρ¨u−div(σ(u)) = βin Ω,(2.2a)

u=uDon ΓD,(2.2b)

σ(u)·n=τon ΓN.(2.2c)

with initial conditions

u(0) = g, ˙u(0) = h.(2.2d)

Therein, udenotes the unknown displacement field with linearized strain tensor ε(u)∈

Rd×d

sym and stress tensor σ(u)∈Rd×d

sym, given by

ε(u) := 1

2∇u+ (∇u)T, σ(u) := λtr ε(u)Id+ 2µε(u).

Note that the stress depends on the material constants λand µ, the so-called Lam´e

parameters, and that Iddenotes the d×didentity matrix whereas tr denotes the trace of

a matrix, i.e., the sum of the diagonal entries. Furthermore, ρ > 0 denotes the constant

density of the material and nthe outer normal vector along the boundary.

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 5

With the inner product for matrices, A:B:= tr(ABT) = Pi,j AijBij, we define the

linear stiffness operator K:V → V∗by

hKu, viV∗,V:= ZΩ

σ(u) : ε(v) dx.(2.3)

By Korn’s inequality [BS08, Ch. 11.2], Kis coercive on VBif ΓDhas positive measure.

Furthermore, Kis symmetric and bounded. Thus, there exist positive constants k1and

k2such that for all u∈ VBand v,w∈ V it holds that

k1kuk2≤ hKu, uiV∗,V,hKv, wiV∗,V≤k2kvkkwk.(2.4)

Note that the symmetry of the operator implies that we may write hKu, uiV∗,V=|K1/2u|2.

2.2.2. Damping. We like to enrich the mathematical model (2.2) by a dissipation term.

Note that the choice of the damping model is a delicate task and depends strongly on the

desired effects. Often viscous damping [Hug87, Ch. 7.2] is considered which corresponds to

a generalization of Hooke’s law. The popular generalization of the mass proportional and

stiffness proportional damping is called Rayleigh damping [CP03, Ch. 12]. Since this quite

common approach has no physical justification [Wil98, Ch. 19], we allow more general

nonlinear damping terms. For this, we define a nonlinear damping operator D:V → V∗

which is assumed to be Lipschitz continuous and strongly monotone, i.e., there exist

constants d0,d1, and d2such that for all u,v∈ V it holds that

kDu− DvkV∗ ≤d2ku−vk, d1ku−vk2−d0|u−v|2≤ hDu− Dv, u −viV∗,V.(2.5)

Furthermore, we may assume w.l.o.g. D(0) = 0, see [ET10a, p. 181], and thus,

kDukV∗≤d2kuk, d1kuk2−d0|u|2≤ hDu, uiV∗,V.

Remark 2.1.Because of the continuous embedding V,→ H, we have | · | ≤ Cembk·k. In

the case d0C2

emb < d1, we can write

hDu, uiV∗,V≥d1kuk2−d0|u|2≥d1−d0C2

embkuk2.

Thus, we may assume either d0= 0 or d0C2

emb ≥d1.

2.2.3. Dirichlet Boundary Conditions. We include the inhomogeneous Dirichlet boundary

conditions in a weak form, i.e., with the help of Lagrange multipliers [Sim00, Sim13].

This leads to a dynamic saddle point problem which is advantageous if the Dirichlet

data depends e.g. on the motion of other bodies as in flexible multibody dynamics.

Furthermore, the considered setting includes boundary control [Tr¨o09].

Within this paper, we denote the trace operator, i.e., the extension of (2.1), by B:V →

Q∗. Note that the space VBequals the kernel of the operator Band let Vcdenote any

complement such that

V=VB⊕ Vc.

Since Vis a Hilbert space, the canonical choice is the orthogonal complement Vc= (VB)⊥V.

In any case, the operator Bis an isomorphism, if restricted to Vc, and Bsatisfies an inf-sup

condition, i.e., there exists a constant β > 0 with

inf

q∈Q

sup

v∈V

hBv, qi

kvk kqkQ

=β > 0,

In other words, the operator Bhas a continuous right inverse which we denote by B−.

The corresponding continuity constant is given by CB−, i.e., kB−· k ≤ CB−k · kQ∗. Finally,

6 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

its dual operator B∗:Q → Vo

B⊆ V∗defines an isomorphism, where Vo

Bdenotes the polar

set (also called annihilator), i.e.,

B:= f∈ V∗| hf, vi= 0 for all v∈ VB.(2.6)

2.3. Formulation as Operator DAE. The weak formulation of equation (2.2a) with

additional damping can be written in operator form. This then equals an operator ODE,

i.e., an ODE in an abstract setting. Including the inhomogeneous boundary conditions

by the Lagrangian method, we add a constraint and thus, we obtain an operator DAE. In

this case, the solution consists of the deformation variable uand the Lagrange multiplier

λ.

We consider two different operator formulations. Either way, we assume for the data of

the right-hand sides F ∈ L2(0, T;V∗) and G ∈ H2(0, T;Q∗). Note that the regularity of G

in the time variable is a necessary condition of the existence of solutions.

2.3.1. Original Formulation. To ensure that the introduced operators are defined for the

solution, we assume that the deformation variable satisfies u∈H1(0, T;V) with second

derivative ¨u∈L2(0, T;V∗). Note that ˙u∈L2(0, T;H) is not sufficient because of the

damping term. As search space for the Lagrange multiplier we consider L2(0, T;Q). Thus,

the dynamic saddle point problem in operator form reads:

Find u∈H1(0, T;V) with ¨u∈L2(0, T;V∗) and λ∈L2(0, T;Q) such that

ρ¨u(t) + D˙u(t) + Ku(t) + B∗λ(t) = F(t) in V∗,(2.7a)

Bu(t) = G(t) in Q∗

(2.7b)

is satisfied for t∈(0, T) a.e. and initial conditions

u(0) = g∈ V,˙u(0) = h∈ H.(2.7c)

Note that the assumed regularity of the solution implies u∈C([0, T]; V) and ˙u∈C([0, T]; H),

cf. [Rou05, Ch. 7]. Thus, the initial conditions are well-posed for g∈ V and h∈ H. As

classical DAEs require consistent initial data, we have to expect a similar condition in the

operator case.

Remark 2.2.Because of the constraint (2.7b), the initial data have to satisfy Bg=G(0).

Thus, we obtain the decomposition g=g0+B−G(0) with g0∈ VBwhich is a consistency

condition for g. For hwe get the decomposition h=h0+B−˙

G(0) with h0∈ H. Note that,

since B−˙

G(0) ∈ Vc,→ H, this does not give a restriction for h.

System (2.7) is called an operator DAE, since a spatial discretization by finite elements

yields a DAE of index 3 [Alt13]. Since high-index DAEs are known to be very sensitive

to perturbations, their numerical approximation is a difficult task and numerical time

integration methods may even diverge [LP86]. For a simulation it is therefore advisable

to perform an index reduction which yields an equivalent system of equations which is of

index one. In the infinite-dimensional case, a similar approach is possible.

2.3.2. Regularized Formulation. The regularization of semi-explicit operator DAEs is pre-

sented in [Alt13, Alt15]. This regularization then results in an equivalent operator DAE

whose spatial discretization leads to a DAE of index 1. Thus, the resulting system is

better suited for numerical integration.

The regularization involves an extension of the system by the so-called hidden con-

straints and additional dummy variables. Furthermore, the deformation variable uis split

into u=u1+u2where u1is the differential part in VBand u2the part of the deformation

which is already fixed due to the given constraint. The system then reads as follows:

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 7

Find u1∈H1(0, T;VB) with ¨u1∈L2(0, T;V∗) as well as u2,v2,w2∈L2(0, T;Vc) and

λ∈L2(0, T;Q) such that

ρ(¨u1+w2) + D( ˙u1+v2) + K(u1+u2) + B∗λ=Fin V∗,(2.8a)

Bu2=Gin Q∗,(2.8b)

Bv2=˙

Gin Q∗,(2.8c)

Bw2=¨

Gin Q∗,(2.8d)

is satisfied for t∈(0, T) a.e. with initial conditions

u1(0) = g0∈ VB,˙u1(0) = h0∈ H.(2.8e)

With the given assumptions on the involved operators from Section 2.2, system (2.8) has

a unique solution (u1, u2, v2, w2, λ), see [Alt13]. Furthermore, it has been shown that the

operator DAE is well-posed in the sense that the map

g0, h0,F,G7→ u1, u2, v2, w2,¨u1+D˙u1+B∗λ

is linear and continuous as mapping

VB× H × L2(0, T;V∗)×H2(0, T;Q∗)→

C([0, T],V)∩C1([0, T],H)×L2(0, T;Vc)3×L2(0, T;V∗).

Remark 2.3.The used regularization technique also applies to flow equations such as the

Stokes or Oseen equations. In [AH13] it has been shown that this formulation is beneficial

for numerical simulations and even allows semi-explicit time integration schemes.

3. Discretization and Stability

To derive a priori error estimates and the convergence proofs, we apply standard tech-

niques from abstract ODE theory as in [ET10a]. For this, we construct piecewise constant

and linear (in time) approximations of the variables of interest. The a priori estimates

then show the boundedness of the approximation independent of the step size such that a

weakly convergent subsequence can be extracted.

3.1. Time Discretization. In the fields of elastodynamics and multibody dynamics, the

Newmark scheme [New59, GC01] as well as further developments like the generalized-α

methods [CH93, AB07] are widely used. However, these schemes are not suitable for

the convergence of operator equations, since also derivatives of the approximations of the

previous time step are used [Eˇ

ST13]. Thus, we restrict ourselves to the scheme which

corresponds to the implicit Euler method applied to the equivalent first-order system.

We only consider equidistant time steps with step size τ. Let ujdenote the approxima-

tion of uat time tj=jτ. For the temporal discretization we then replace the derivatives

˙uand ¨uat time tjby

˙u(tj)→uj−uj−1

τ=: Duj,¨u(tj)→uj−2uj−1+uj−2

τ2=: D2uj.

The convergence of this scheme for index-3 DAEs, i.e., for the finite-dimensional setting

arising in multibody dynamics, is discussed in [LP86]. We emphasize that the analysis

used in [LP86] assumes that the constraint is solved with high accuracy, namely up to the

order of O(τ3). This is not necessary if the index of the system is reduced first.

8 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

3.1.1. Function Evaluations. For the (formal) application of a discretization scheme to an

operator equation, we need function evaluations of the right-hand sides. However, the

given data may not be continuous. Thus, function evaluations of the right-hand side have

to be replaced, e.g., by an integral mean over one time step.

Consider a Bochner integrable function F ∈ L2(0, T ;X) with a real Banach space X

and an equidistant partition 0 = t0< t1<· · · < tn=Tof [0, T ]. We define Fj∈Xby

the Bochner integral over one time step τ, i.e.,

Fj:= 1

τZtj

tj−1

F(t) dt∈X.

Note that this is well-defined for F ∈ L2(0, T;X). In this way, we define the piecewise

constant (abstract) function Fτ: [0, T]→Xby

Fτ(t) := Fjfor t∈(tj−1, tj](3.1)

and a continuous extension in t= 0. An easy calculation shows that Fτ∈L2(0, T;X)

satisfies the inequality

kFτk2

L2(0,T;X)=τ

j=1

kFjk2

X≤

j=1 Ztj

tj−1

kF(t)k2

Xdt=kFk2

L2(0,T;X).(3.2)

One important property of Fτis the strong convergence to F.

Lemma 3.1 ([Tem77, Ch. III, Lem. 4.9]).Consider F ∈ L2(0, T;X)with its approxima-

tion Fτas defined in (3.1). Then, Fτ→ F in L2(0, T;X), i.e., kFτ− FkL2(0,T;X)→0,

as τ→0.

For continuous functions F ∈ C([0, T]; X) function evaluations are well-defined. In this

case, we may define

Fτ(t) := F(tj)∈Xfor t∈(tj−1, tj].(3.3)

Again we consider a continuous extension in t= 0 and obtain Fτ→ F in L2(0, T ;X).

If F ∈ H1(0, T;X), then we discretize Fby means of function evaluations as in (3.3)

and ˙

Fby the integral mean as in (3.1). This approach has the nice property that the

discrete derivative of Fjequals the approximation of the derivative, i.e.,

DFj=Fj− Fj−1

τ=˙

Fj.

3.1.2. Semi-discrete Equations. Replacing the derivatives by discrete derivatives, i.e., ˙u(tj)

by Dujand ¨u(tj) by D2uj, we obtain from the differential equation

ρD2uj

1+wj

2+DDuj

1+vj

2+Kuj

1+uj

2+B∗λj=Fj.(3.4a)

This equation has to be solved for j= 2, . . . , n and is still stated in the dual space of

Vand thus, equals a PDE in the weak formulation. The three constraints (2.8b)-(2.8d)

result in

Buj

2=Gj,Bvj

2=˙

Gj,Bwj

2=¨

Gjin Q∗.(3.4b)

Remark 3.1 (Special case G ≡ 0).Consider the case where Gvanishes, i.e., the homo-

geneous Dirichlet case. This directly implies uj

2=vj

2=wj

2= 0 and thus, the problem

reduces to an operator ODE on the kernel VB, namely

ρD2uj

1+DDuj

1+Kuj

1=Fjin V∗

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 9

Before we derive stability results for the discrete approximations, we have to discuss

the solvability of the semi-discrete system (3.4).

Lemma 3.2. With the assumptions introduced in Section 2, system (3.4) has a unique

solution (uj

1, uj

2, vj

2, wj

2, λj)for each time step j= 2, . . . , n if the step size satisfies τ < ρ/d0.

In the case d0= 0, there is no step size restriction.

Proof. The invertibility of the operator Bin Vcimplies that the three equations in (3.4b)

give unique approximations uj

2,vj

2, and wj

2, respectively. Consider equation (3.4a) re-

stricted to test functions in VB. We define the operator A:VB→ V∗

Band the functional

Fj∈ V∗by

Au:= ρ

τ2u+Du−uj−1

τ+vj

2+Ku, ¯

Fj:= Fj+ρ

τ22uj−1

1−uj−2

1−ρwj

2− Kuj

Then, equation (3.4a) can be written in the form Auj

1=¯

Fjin V∗

B. Obviously, the operator

Ais continuous. Using (2.4) and (2.5), we also have that Ais monotone, since

Au− Av, u −v≥ρ

τ2|u−v|2+d1

τku−vk2−d0

τ|u−v|2+k1ku−vk2

=d1/τ +k1ku−vk2+ρ/τ2−d0/τ|u−v|2.

This shows that hAu− Av, u −vi ≥ k1ku−vk2for τ < ρ/d0and thus, the existence

of a solution uj

1∈ VBusing the Browder-Minty theorem [GGZ74, Ch. III, Th. 2.1]. The

strong monotonicity of Aalso implies the uniqueness of the solution. Finally, the unique

solvability for λjfollows from the invertibility of B∗:Q→Vo

B, cf. Section 2.2.3. 

3.2. Stability Estimates. Within this subsection, we use the abbreviation

1:= Duj

1=uj

1−uj−1

τ.

Furthermore, we assume u1

1and v1

1to be the fixed initial data of the semi-discrete problem

(3.4), i.e., approximations of the initial data u1(0) = g0and ˙u1(0) = h0. Note that this

also defines u0

1which - in the limit - coincides with u1

In the following lemma we give a stability estimate of the semi-discrete approximations.

Note that this includes a step size restriction due to the nonlinear damping term.

Lemma 3.3 (Stability).Assume right-hand sides F ∈ L2(0, T;V∗),G ∈ H2(0, T;Q∗)and

initial approximations u1

1∈ VB,v1

1∈ H. Let the approximations uj

1,uj

2,vj

2, and wj

2be

given by the semi-discrete system (3.4) and let the step size satisfy τ < ρ/8d0. Then, there

exists a constant c > 0such that for all k≥2the inequality

ρvk

1

2+ρ

j=2 vj

1−vj−1

1

2+τd1

j=2 

vj

1



2+k1

uk

1



2≤c28d0T/ρM2

(3.5)

is satisfied with a constant M=hku1

1k2+|v1

1|2+kFk2

L2(0,T;V∗)+kGk2

H2(0,T;Q∗)i1/2.

Proof. The equations in (3.4b) directly lead to the estimates

kuj

2k ≤ CB−kGjkQ∗,kvj

2k ≤ CB−k˙

GjkQ∗,kwj

2k ≤ CB−k¨

GjkQ∗.(3.6)

The remainder of the proof follows the ideas of the proof of [ET10a, Th. 1] although a

different time discretization scheme is used. We only consider the case d0>0, since the

10 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

proof for d0= 0 works in the same manner but with less difficulties. Within the proof, we

take several times advantage of the equality

2(a−b)a=a2−b2+ (a−b)2.(3.7)

We test equation (3.4a) with the discrete derivative vj

1∈ VB,j≥2. This leads to

ρDvj

1, vj

1+Dvj

1+vj

2, vj

1+Kuj

1, vj

1=Fj, vj

1−ρwj

2, vj

1−Kuj

2, vj

1.(3.8)

For the terms on the left-hand side, we estimate separately

ρDvj

1, vj

1=ρ

τvj

1−vj−1

1, vj

1(3.7)

=ρ

2τhvj

1

2−vj−1

1

2+vj

1−vj−1

1

2i,

for the damping term

Dvj

1+vj

2, vj

1=Dvj

1+vj

2− Dvj

2, vj

1+Dvj

2, vj

1

(2.5)

≥d1kvj

1k2−d0|vj

1|2−d2kvj

1kkvj

≥d1kvj

1k2−d0|vj

1|2−d1

6kvj

1k2−3d2

2d1

kvj

2k2,

and finally, for the stiffness term

Kuj

1, vj

1=1

τKuj

1, uj

1−uj−1

1(3.7)

≥1

2τK1/2uj

1

2−1

2τK1/2uj−1

1

Using the Cauchy-Schwarz inequality, followed by an application of Youngs inequality

[Eva98, App. B], for the right-hand side of (3.8) we obtain

Fj, vj

1−ρwj

2, vj

1−Kuj

2, vj

1

≤ kFjkV∗kvj

1k+ρ|wj

2||vj

1|+k2kuj

2kkvj

≤3

2d1

kFjk2

V∗+d1

6kvj

1k2+ρ2

4d0

|wj

2|2+d0|vj

1|2+3k2

2d1

kuj

2k2+d1

6kvj

1k2.

Thus, a multiplication of (3.8) by 2τimplies with the estimates above that

ρh|vj

1|2− |vj−1

1|2+|vj

1−vj−1

1|2i+τd1kvj

1k2−4τd0|vj

1|2+K1/2uj

1

2−K1/2uj−1

1

≤τ"3

kFjk2

∗+3k2

kuj

2k2+3d2

kvj

2k2+ρ2

2d0

|wj

2|2#.(3.9)

With the estimates of uj

2,vj

2, and wj

2from equation (3.6) we can bound the right-hand

side of (3.9) by cτkFjk2

V∗+kGjk2

Q∗+k˙

Gjk2

Q∗+k¨

Gjk2

Q∗. Here, c > 0 denotes a generic

constant which depends on CB−,ρ,d0,d1,d2, and k2.

Before we sum over jand make benefit of several telescope sums, we have to deal with

the term 4τd0|vj

1|2on the left-hand side of (3.9). For this, we use arguments which are

used to prove discrete versions of the Gronwall lemma [Emm99]. With κ:= 4d0/ρ and

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 11

aj:= (1 −κτ)j, we estimate

ρhaj|vj

1|2−aj−1|vj−1

1|2+aj−1|vj

1−vj−1

1|2i+τd1aj−1kvj

1k2+ajK1/2uj

1

2−aj−1K1/2uj−1

1

=aj−1hρ(1 −κτ)|vj

1|2−ρ|vj−1

1|2+ρ|vj

1−vj−1

1|2+τd1kvj

1k2

+ (1 −κτ)K1/2uj

1

2−K1/2uj−1

1

≤aj−1hρ|vj

1|2−ρ|vj−1

1|2+ρ|vj

1−vj−1

1|2+τd1kvj

1k2−4τd0|vj

1|2

+K1/2uj

1

2−K1/2uj−1

1

(3.9)

≤aj−1τc kFjk2

V∗+kGjk2

Q∗+k˙

Gjk2

Q∗+k¨

Gjk2

Q∗.

Note that we have used the fact that, due to the assumption on the step size τ, we have

0< aj<1 for all j≥1 and κ≥0. The summation of this estimate for j= 2, . . . , k then

yields

ρakvk

1

2+ρ

j=2

aj−1vj

1−vj−1

1

2+τd1

j=2

aj−1kvj

1k2+akK1/2uk

1

≤ρa1v1

1

2+a1K1/2u1

1

2+τc

j=2

aj−1kFjk2

V∗+kGjk2

Q∗+k˙

Gjk2

Q∗+k¨

Gjk2

Q∗.

Finally, we divide by akand use the estimates aj> akfor j < k and a−k≤4κT . The

latter inequality follows from 1/τ > 2κand the monotonicity of the sequence (1 + x/n)n

ak= (1 −κτ)k= (1 −κT/n)k>(1 −κT/n)n≥(1 −1/2)2κT = 4−κT .

This then leads to the final result

ρvk

1

2+ρ

j=2 vj

1−vj−1

1

2+τd1

j=2 

vj

1



2+k1

uk

1



≤4κT nρv1

1

2+k2

u1

1



2+τc

j=2 kFjk2

V∗+kGjk2

Q∗+k˙

Gjk2

Q∗+k¨

Gjk2

Q∗o.

With the stability estimate (3.5) in hand, we are able to show the uniform boundedness

of the approximation sequences defined in the following subsection.

3.3. Global Approximations. In this subsection, we define global approximations of

u1,u2,v2, and w2. First, we define U1,τ ,ˆ

U1,τ : [0, T]→ VBby

U1,τ (t) := uj

1,ˆ

U1,τ (t) := uj

1+ (t−tj)vj

for t∈(tj−1, tj] and j≥2 with U1,τ ≡ˆ

U1,τ ≡u1

1on [0, t1]. By the stability estimate (3.5)

of Lemma 3.3 we directly obtain the uniform boundedness of U1,τ and ˆ

U1,τ in L∞(0, T;VB).

Thus, there exists a weak limit U1∈L∞(0, T;VB) with U1,τ ,ˆ

U1,τ

∗

−* U1in L∞(0, T;VB) as

well as U1,τ ,ˆ

U1,τ * U1in L2(0, T;VB). Note that the limits of the two sequences coincide,

since



ˆ

U1,τ −U1,τ 



L2(0,T;H)=

j=1 Ztj

tj−1(t−tj)vj

1

2dt≤

j=1

τ3vj

1

2(3.5)

≤cτ2M2→0.

12 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

From this, the agreement of the limits in L2(0, T;VB) follows by the assumed embedding

VB,→ H given by the Gelfand triple.

In an analogous way, we define the piecewise constant functions U2,τ ,V2,τ ,W2,τ : [0, T]→

Vc. We set

U2,τ (t) := uj

2, V2,τ (t) := vj

2, W2,τ (t) := wj

for t∈(tj−1, tj] and j≥1 with a continuous extension in t= 0. By equation (3.4b) we

have BU2,τ =Gτ,BV2,τ =˙

Gτ, and BW2,τ =¨

Gτ. Thus, Lemma 3.1 implies that

U2,τ →U2, V2,τ →V2, W2,τ →W2in L2(0, T, V).

Note that the limits U2,V2, and W2solve the equations BU2=G,BV2=˙

G, and BW2=¨

respectively. This means nothing else than the (strong) convergence of U2,τ ,V2,τ , and W2,τ

to the solutions of (2.8b)-(2.8d).

Finally, we define two different approximations of the velocity in form of a piecewise

constant and a piecewise linear approximation, namely

V1,τ (t) := vj

1,ˆ

V1,τ (t) := vj

1+ (t−tj)Dvj

for t∈(tj−1, tj] and j≥2 with V1,τ ≡ˆ

V1,τ ≡v1

1on [0, t1]. An illustration is given in

Figure 3.1.

t1t2t3

V1,τ

Figure 3.1. Illustration of the global approximations V1,τ and ˆ

V1,τ of ˙u1.

For the piecewise constant approximation, we obtain the estimate

kV1,τ k2

L2(0,T;V)=ZT

kV1,τ (t)k2dt=τ

j=1 

vj

1



2(3.5)

≤τkv1

1k2+cM2.

Up to now, we have only assumed v1

1∈ H. In order to obtain a uniform bound of V1,τ , we

have to assume v1

1∈ V. This then implies the existence of a weak limit V1∈L2(0, T;VB),

i.e., V1,τ * V1in L2(0, T;VB). In the same manner we obtain a bound of the piecewise

linear approximation,



ˆ

V1,τ 



L2(0,T;V)=τkv1

1k2+

j=2 Ztj

tj−1

vj

1+ (t−tj)Dvj

1



2dt≤4τ

j=1 

vj

1



As before, we show that V1,τ and ˆ

V1,τ have the same limit V1. For this, by Lemma 3.3 we

calculate that



ˆ

V1,τ −V1,τ 



L2(0,T;H)=

j=1 Ztj

tj−1ˆ

V1,τ (t)−V1,τ (t)

2dt≤τ

j=2 vj

1−vj−1

1

2≤τcM2→0.

In the following, we show that the limit function V1equals the derivative of U1in the

generalized sense. For this, we use the limits ˆ

U1,τ * U1and V1,τ * V1in L2(0, T;VB).

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 13

Note, that d

dtˆ

U1,τ =V1,τ a.e. but not in the interval [0, τ]. Applying the integration by

parts formula with an arbitrary functional f∈ V∗

Band Φ ∈C∞

0([0, T]), we see that

0f, U1˙

Φ dt= lim

τ→0ZT

0f, ˆ

U1,τ ˙

Φ dt=−lim

τ→0ZT

0f, ˙

U1,τ Φ dt

=−lim

τ→0ZT

0f, V1,τ Φ dt−Zτ

0f, v1

1Φ dt=−ZT

0f, V1Φ dt.

Note that the integral over [0, τ] vanishes in the limit, since the integrand is bounded in-

dependently of the step size. As a result, the limit function U1has a generalized derivative

and ˙

U1=V1∈L2(0, T;VB).

Finally, we mention that also D(V1,τ +V2,τ ) gives a uniformly bounded sequence in

L2(0, T;V∗) due to the continuity of the damping operator D. Thus, there exists a weak

limit D∈L2(0, T;V∗) with

D(V1,τ +V2,τ )*Din L2(0, T;V∗).

One aim of the next subsection is to show that aequals D(V1+V2), i.e., the limit of the

damping term equals the damping operator applied to the limit functions.

4. Convergence

This section is devoted to the analysis of the limiting behaviour of the discrete approx-

imations. Furthermore, we analyse the influence of additional perturbations which then

shows the convergence of the Rothe method applied to the operator DAE (2.8).

4.1. Deformation Variable. In order to pass to the limit with τ→0 it is beneficial

to rewrite equation (3.4a) in terms of the global approximations. In this subsection, we

only consider test functions in VBin order to eliminate the Lagrange multiplier from the

system. The semi-discrete system has the form

ρ˙

V1,τ +W2,τ +DV1,τ +V2,τ +KU1,τ +U2,τ =Fτin V∗

(4.1)

for t∈(τ, T ) a.e.. Writing equation (4.1) in its actual meaning with test functions v∈ VB,

Φ∈C∞

0([0, T]) and applying the integration by parts formula once, we get

−ρˆ

V1,τ , v˙

Φ + ρW2,τ , vΦ + DV1,τ +V2,τ , vΦ + KU1,τ +U2,τ , vΦ dt

=ZT

0Fτ, vΦ dt.

Passing to the limit, by the achievements of the previous section we obtain that

0ρV1, v˙

Φ dt=ZT

0ρW2+D+KU1+U2− F, vΦ dt.

Recall that Ddenotes the weak limit of D(V1,τ +V2,τ ) in L2(0, T;V∗). This implies that

V1has a generalized derivative ˙

V1∈L2(0, T;V∗

B) which satisfies the equation

ρ˙

V1+ρW2+D+K(U1+U2) = Fin V∗

B.(4.2)

The remaining part of this subsection is devoted to the proof that the weak limits U1,U2,

V2, and W2solve the operator DAE (2.8a) in V∗

B. With equation (4.2) at hand, it remains

to show that Dequals D(V1+V2). In order to show this, we present two preparatory

lemmata.

14 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

Lemma 4.1. For t=Tthe sequence ˆ

V1,τ satisfies ˆ

V1,τ (T)* V1(T)in H. Furthermore,

we obtain the estimate lim infτ→0˙

V1,τ , V1,τ ≥˙

V1, V1.

Proof. We follow the idea of the proof of [ET10b, Th. 5.1], adapted to the given operator

equation. First we show that ˆ

V1,τ (T)* V1(T) in Has well as ˆ

V1,τ (0) = V1(0). Because

of the stability estimate in Lemma 3.3, the final approximation ˆ

V1,τ (T) = vn

1is uniformly

bounded in H. Thus, there exists a weak limit ξ∈ H which satisfies

1=ˆ

V1,τ (T)* ξ in H.

Through the integration by parts formula and with w∈ VBand Φ ∈C1([0, T]), we obtain

ρV1(T), wΦ(T)−ρV1(0), wΦ(0)

=ρ˙

V1, wΦ+ρV1, w ˙

Φ

(4.2)

=F − ρW2−D− K(U1+U2), wΦ+ρV1, w ˙

Φ

(4.1)

=F − Fτ, wΦ−ρW2−W2,τ , wΦ−D− D(V1,τ +V2,τ ), wΦ

−K(U1+U2)− K(U1,τ +U2,τ ), wΦ+ρV1, w ˙

Φ+ρ˙

V1,τ , wΦ

=F − Fτ, wΦ−ρW2−W2,τ , wΦ−D− D(V1,τ +V2,τ ), wΦ

−K(U1+U2)− K(U1,τ +U2,τ ), wΦ+ρV1−ˆ

V1,τ , w ˙

Φ

+ρˆ

V1,τ (T), wΦ(T)−ρˆ

V1,τ (0), wΦ(0)

→ρξ, wΦ(T)−ρv1

1, wΦ(0).

Thus, we have vn

1=ˆ

V1,τ (T)* ξ =V1(T) in Hand V(0) = v1

1. Note that at this point we

need that the embedding VB,→ H is dense. A direct consequence of the weak convergence

is that |V1(T)| ≤ lim infτ→0|vn

1|. With the calculation

˙

V1,τ , V1,τ =

j=1 vj

1−vj−1

1, vj

1≥ −1

j=1 vj−1

1

2−vj

1

2=1

2|vn

1|2−1

2|v1

1|2

we finally conclude

lim inf

τ→0˙

V1,τ , V1,τ ≥1

2lim inf

τ→0|vn

1|2− |v1

1|2≥1

2V1(T)

2−1

2V1(0)

2=˙

V1, V1.

Remark 4.1.The fact that ˆ

V1,τ (T)* V1(T) in Hand ˆ

V1,τ (0) = V1(0), as shown in

Lemma 4.1, implies that for w∈ VBand Φ ∈C2([0, T]) it holds that

lim

τ→0˙

V1,τ , w ˙

Φ= lim

τ→0−ˆ

V1,τ , w¨

Φ+ˆ

V1,τ (T), w˙

Φ(T)−ˆ

V1,τ (0), w˙

Φ(0)

=−V1, w¨

Φ+V1(T), w˙

Φ(T)−V1(0), w˙

Φ(0) = ˙

V1, w ˙

Φ.

In the following lemma, we consider the stiffness operator K.

Lemma 4.2. The sequences U1,τ ,U2,τ , and V1,τ satisfy the estimate

lim inf

τ→0K(U1,τ +U2,τ ), V1,τ ≥K(U1+U2), V1.

Proof. Because of the linearity of Kand the strong convergence of U2,τ it is sufficient to

analyse the lim inf of hKU1,τ , V1,τ iand show that lim infτ→0hKU1,τ , V1,τ i ≥ hKU1, V1i. For

this, we proceed as in the proof of Lemma 4.1.

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 15

Lemma 3.3 implies the boundedness of ˆ

U1,τ (T) = un

1in Vsuch that there exists an

element ξ∈ VBwith ˆ

U1,τ (T)* ξ in VB. We show that K1/2ξ=K1/2U1(T) and K1/2u1

K1/2U1(0). Using the limit equation (4.2) and the semi-discrete equation (4.1) with test

functions w∈ VBand Φ ∈C2([0, T]), we obtain

KU1(T),wΦ(T)−KU1(0), wΦ(0)

=K˙

U1, wΦ+KU1, w ˙

Φ

(4.2)

=K˙

U1, wΦ+F − ρW2−D− KU2−ρ˙

V1, w ˙

Φ

(4.1)

=F − Fτ, w ˙

Φ−ρW2−W2,τ , w ˙

Φ−D− D(V1,τ +V2,τ ), w ˙

Φ

−KU2− KU2,τ , w ˙

Φ−ρ˙

V1−˙

V1,τ , w ˙

Φ+K˙

U1, wΦ+KU1,τ , w ˙

Φ.

Passing to the limit, we make use of Remark 4.1 which implies that the term including

V1vanishes. In addition, we use the fact that, passing to the limit, we may replace U1,τ

by ˆ

U1,τ , since they have the same limit. Thus, another application of the integration by

parts formula then leads to

KU1(T), wΦ(T)−KU1(0), wΦ(0) = Kξ, wΦ(T)−Ku1

1, wΦ(0).

Since hK·,·i defines an inner product in VB, we conclude that U1(T) = ξand U1(0) = u1

1in

VB. As a result, we obtain K1/2un

1*K1/2ξ=K1/2U1(T) in Hand K1/2u1

1=K1/2U1(0).

Since U1,τ and V1,τ are both piecewise linear, as in the proof of Lemma 4.1, we calculate

that

KU1,τ , V1,τ =

j=1 Kuj

1, uj

1−uj−1

1≥1

2Kun

1, un

1−1

2Ku1

1, u1

1+τKu1

1, v1

1

2K1/2un

1

2−1

2K1/2u0

1

2+τKu1

1, v1

1.

Note that the term τhKu1

1, v1

1ivanishes as τ→0, since u1

1and v1

1are fixed. By the

property |K1/2U1(T)| ≤ lim infτ→0|K1/2un

1|we finally summarize the partial results to

lim inf

τ→0KU1,τ , V1,τ ≥lim inf

τ→0

2K1/2un

1

2−1

2K1/2u1

1

≥1

2K1/2U1(T)

2−1

2K1/2U1(0)

2=KU1,˙

U1=KU1, V1.

With the previous two lemmata we are now able to prove that the limit of the damping

term equals the damping operator applied to the limit functions.

Theorem 4.3. Consider problem (2.8) with right-hand sides F ∈ L2(0, T;V∗),G ∈

H2(0, T;Q∗)and initial approximations u1

1=g0,v1

1=h0∈ VB. Then, we have D=

D(V1+V2)and the (weak) limits U1,U2,V2, and W2solve the operator equation (2.8a)

for test functions v∈ VB.

Proof. We consider the semi-discrete equation (4.1) tested by V1,τ and subtract the term

hD(V1,τ +V2,τ )− Dw, V1,τ +V2,τ −wiwith w∈L2(0, T;V), which is non-negative because

of the monotonicity of the damping operator, cf. Section 2.2.2. This then leads to

0≥ρ˙

V1,τ , V1,τ +ρW2,τ , V1,τ +KU1,τ +U2,τ , V1,τ −Fτ, V1,τ 

+DV1,τ +V2,τ , w −V2,τ +Dw, V1,τ +V2,τ −w.

16 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

The application of the lim inf on both sides in combination with Lemmata 4.1 and 4.2

then leads to

0≥ρ˙

V1, V1+ρW2, V1+K(U1+U2), V1−F, V1

+D, w −V2+Dw, V1+V2−w.

Note that we have used the fact that the sequences V2,τ and W2,τ converge strongly in

L2(0, T;V) and that Dequals the weak limit of D(V1,τ +V2,τ ). Rearranging the terms and

applying the limit equation (4.2), we then obtain

Dw, w −V1−V2≥ρ˙

V1+ρW2+K(U1+U2)− F, V1+D, w −V2

(4.2)

=−D, V1+D, w −V2

=D, w −V1−V2.

Following the Minty trick [RR04, Lem. 10.47], i.e., choosing w:= V1+V2±sv with an

arbitrary function v∈L2(0, T;V) and s∈[0,1], we conclude that D=D(V1+V2). Thus,

with V1=˙

U1the limit equation (4.2) turns to

ρ¨

U1+ρW2+D(˙

U1+V2) + K(U1+U2) = Fin V∗

It remains to check whether U1satisfies the initial conditions. Note that ˙

U1(0) = V1(0) =

1=h0was shown within the proof of Lemma 4.2, whereas U1(0) = u1

1=g0was proved

in Lemma 4.1. 

4.2. Lagrange Multiplier. Up to this point, the obtained convergence results exclude

the Lagrange multiplier, since we only considered test functions in the kernel of the con-

straint operator B. To analyse the limiting behaviour of the Lagrange multiplier, we test

equation (3.4a) by functions v∈ Vc. In terms of the global approximations and with

Λτ(t) := λjfor t∈(tj−1, tj], this equation can be written in the form

ρ˙

V1,τ +W2,τ +DV1,τ +V2,τ +KU1,τ +U2,τ +B∗Λτ=Fτin V∗.(4.3)

Unfortunately, the given setting does not allow to find a uniform bound for Λτin L2(0, T;Q).

The reason is the absence of an upper bound of τPn

j=1 kDvj

1k2

V∗. We obtain this bound

only within the weaker norm of V∗

B. However, we show that the primitive of Λτ, namely

Λτ(t) := Rt

0Λ(s) ds, converges to the solution of the considered operator DAE in a weaker

sense.

In order to obtain an equation for ˜

Λτ, we have to integrate equation (4.3) over the

interval [0, t]. For an arbitrary test function v∈ V, this then leads to the equation

ρˆ

V1,τ +˜

W2,τ , v+˜

D, v+K˜

U1,τ +˜

U2,τ , v+B∗˜

Λτ, v=˜

Fτ, v+ρv1

1, v.

Therein, ˜

Fτ,˜

U1,τ ,˜

U2,τ , and ˜

W2,τ denote the integrals of Fτ,U1,τ ,U2,τ , and W2,τ , respec-

tively, and

˜

D(t), v:= Zt

0DV1,τ (s) + V2,τ (s), vds.

Note that the term ρv1

1=ρˆ

V1,τ (0) occurs due to the integration of ˙

V1,τ .

We show that ˜

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

L2(0, T;V∗) which implies that its primitive ˜

Fτis uniformly bounded in C([0, T]; V∗).

Furthermore, we have shown the boundedness of U1,τ ,U2,τ , and W2,τ in L2(0, T;V) in

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 17

Section 3.3. Thus, their primitives ˜

U1,τ ,˜

U2,τ , and ˜

W2,τ are bounded in C([0, T]; V). With

the Cauchy-Schwarz inequality, we calculate

max

t∈[0,T]˜

D(t), v

(2.5)

≤d2ZT

kV1,τ (s) + V2,τ (s)kkvkds≤d2T1/2kV1,τ +V2,τ kL2(0,T;V)kvk.

Recall that the boundedness of V1,τ +V2,τ in L2(0, T;V) was already shown in Section 3.3.

Finally, the estimate

max

t∈[0,T]|ˆ

V1,τ (t)| ≤ max

j|vj

1|(3.5)

≤c24d0T/ρM

shows



˜

Λτ

C([0,T];Q)≤1

βmax

t∈[0,T]

sup

v∈V B∗˜

Λτ(t), v

kvk,

where βis the inf-sup constant. As a result, there exists a limit function ˜

Λ∈Lp(0, T;Q)

such that ˜

Λτ*˜

Λ in Lp(0, T;Q) for all 1 <p<∞. This then leads to the following

convergence result.

Theorem 4.4. Consider problem (2.8) with right-hand sides F ∈ L2(0, T;V∗),G ∈

H2(0, T;Q∗)and initial data u1

1=g0,v1

1=h0∈ VB. Then, the weak limit ˜

Λ∈L2(0, T;Q)

of the sequence ˜

Λτand U1,U2,V2, and W2solve the operator DAE (2.8) in the weak

distributional sense, meaning that for all v∈ V and Φ∈C∞

0([0, T]) it holds that

−ρ˙

U1, v˙

Φ + ρW2+D˙

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

Λ, v˙

Φdt= 0

as well as BU2=G,BV2=˙

G, and BW2=¨

Gin Q∗. Furthermore, U1satisfies the initial

conditions U1(0) = g0and ˙

U1(0) = h0.

Proof. Considering once more equation (4.3) and integrating by parts, for all v∈ V and

Φ∈C∞

0([0, T]) we obtain

−ρˆ

V1,τ , v˙

Φ+ρW2,τ +DV1,τ +V2,τ +KU1,τ +U2,τ −Fτ, vΦ−B∗˜

Λτ, v˙

Φ dt= 0

By the weak convergence of ˜

Λτand the linearity of B∗, we conclude that

0B∗˜

Λτ, v˙

Φ dt→ZT

0B∗˜

Λ, v˙

Φ dt.

The convergence of the remaining terms - also for test functions v∈ V - as well as the

satisfaction of the initial conditions was already shown in Theorem 4.3. 

In summary, we have shown the strong convergence for u2,v2, and w2, the weak conver-

gence for the differential variable u1, and the convergence in the weak distributional sense

for the Lagrange multiplier λ. This result emphasizes that the Lagrange multiplier behaves

qualitatively different than the deformation variables. This is no surprise since already in

the finite-dimensional DAE case one can observe a different behaviour of differential and

algebraic variables [Arn98].

18 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

4.3. Influence of Perturbations. To show the convergence of the deformation variables

and the Lagrange multiplier, we have always assumed that the stationary PDEs were solved

exactly in every time step. Thinking of the Rothe method, where the solution of these

PDEs would only be approximated, e.g. by the finite element method, additional errors

have to be taken into account. Because of this, we analyse in this subsection the influence

of perturbations in the right-hand sides.

We consider perturbations δj∈ V∗as well as θj,ξj,ϑj∈ Q∗, i.e., we solve system (3.4)

with right-hand sides Fj+δj,Gj+θj,˙

Gj+ξj, and ¨

Gj+ϑjinstead of Fj,Gj,˙

Gj, and ¨

Gj.

We denote the solution of the perturbed problem by (ˆuj

1,ˆuj

2,ˆvj

2,ˆwj

2,ˆ

λj). The differences of

the exact and perturbed solution are then given by

1:= ˆuj

1−uj

1, ej

2:= ˆuj

2−uj

2, ej

v:= ˆvj

2−vj

2, ej

w:= ˆwj

2−wj

2.(4.4)

The initial errors in u1

1and v1

1are denoted by e1

1and ˙e1

1, respectively.

Remark 4.2.In some cases, the spatial error of a finite element discretization can be

viewed as a perturbation of the semi-discrete system. Note that the results of this sec-

tion only apply if ej

1∈ VB, i.e., if we consider conforming methods. If this is the case,

then the residuals may be interpreted as perturbations of the right-hand sides, cf. [Alt15,

Sect. 10.4.2].

Considering test functions in VB, the errors in (4.4) satisfy the equation

ρD2ej

1+ρej

w+Dˆvj

1+ ˆvj

2− Dvj

1+vj

2+Kej

1+ej

2=δj.(4.5)

Furthermore, the errors ej

2,ej

v, and ej

wsatisfy the equations Bej

2=θj,Bej

v=ξj, and

Bej

w=ϑjin Q∗which directly yields

kej

2k ≤ CB−kθjkQ∗,kej

vk ≤ CB−kξjkQ∗,kej

wk ≤ CB−kϑjkQ∗.

From equation (4.5) we deduce an estimate of the resulting error ej

1. For this, we follow

again the lines of Lemma 3.3 and test the equation with Dej

1. The only difference takes

place in the estimate of the damping term for which we obtain

Dˆvj

1+ˆvj

2− Dvj

1+vj

2, Dej

1

=Dˆvj

1+ ˆvj

2− Dvj

1+vj

2, Dej

1+ej

v−Dˆvj

1+ ˆvj

2− Dvj

1+vj

2, ej

v

(2.5)

≥d1

Dej

1+ej

v



2−d0Dej

1+ej

v

2−d2

Dej

1+ej

v



ej

v

.

Following the remaining parts of the proof of Lemma 3.3, for k≥2 and sufficiently small

step size τwe yield an estimate of the form

ρDek

1

2+ρ

j=2 Dej

1−Dej−1

1

2+τd1

j=2 

Dej

1+ej

v



2+k1

ek

1



2≤c28d0T/ρM2

The constant Methen includes the initial errors as well as the perturbations. More

precisely, assuming perturbations of comparable magnitude, i.e., δj≈δ,θj≈θ,ξj≈ξ,

and ϑj≈ϑ, we have

e=ke1

1k2+|˙e1

1|2+Thkδk2

V∗

B+kθk2

Q∗+kξk2

Q∗+kϑk2

Q∗i.(4.6)

Summarizing the estimates of this subsection, we obtain the following theorem.

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 19

Theorem 4.5. Consider system (3.4) and the assumptions of Lemma 3.3. Furthermore,

let the perturbations δj∈ V∗and θj,ξj,ϑj∈ Q∗be of the same order of magnitude. Then,

with the constant Mefrom (4.6) and a sufficiently small step size τ, the errors ek

1,ek

2,ek

and ek

wsatisfy

kek

1k2+kek

2k2+kek

vk2+kek

wk2≤ce4d0T/ρM2

This theorem shows that the errors due to perturbations of the right-hand sides are

bounded by these perturbations. Note that this is only true for the regularized operator

DAE (2.8). If we consider the original formulation (2.7) instead, then the error ek

1gets

amplified by a factor 1/τ2, since ξjhas to be replaced roughly by τ−1θjand ϑjby τ−2θj.

Note furthermore that Theorem 4.5 does not include the error in the Lagrange multi-

plier. As already seen in the previous subsection, we are not able to find bounds for the

Lagrange multiplier in the given setting. In the linear case, however, this is possible if we

assume more regularity of the perturbations such as δ∈ H∗, cf. [Alt15].

5. Conclusion

We have shown that the Rothe method, which is very popular in the finite element

community for solving time-dependent PDEs, can also be applied to (regularized) opera-

tor DAEs, i.e., if we include additional constraints. Similar to the finite-dimensional case,

where it is advisable to consider index-1 formulations, we have used the regularized formu-

lation of the operator DAE. With a splitting of the deformation variable into a differential

part and a constrained part, we were able to use PDE techniques to prove the convergence

of the Euler scheme.

Ongoing research also considers higher-order Runge-Kutta methods. We hope that in

this case, under sufficient smoothness assumptions, also the convergence of the Lagrange

multiplier can be shown.

References

[AB07] M. Arnold and O. Br¨uls. Convergence of the generalized-αscheme for con-

strained mechanical systems. Multibody Syst. Dyn., 18(2):185–202, 2007.

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

edition, 2003.

[AH13] R. Altmann and J. Heiland. Finite element decomposition and minimal ex-

tension for flow equations. Preprint 2013–11, Technische Universit¨at Berlin,

Germany, 2013. accepted for publication in M2AN.

[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.

[Arn98] M. Arnold. Zur Theorie und zur numerischen L¨osung von Anfangswertprob-

lemen f¨ur differentiell-algebraische Systeme von h¨oherem Index. VDI Verlag,

D¨usseldorf, 1998.

[Bau10] O. A. Bauchau. Flexible Multibody Dynamics. Solid Mechanics and Its Applica-

tions. Springer-Verlag, 2010.

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

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

[CH93] J. Chung and G. M. Hulbert. A time integration algorithm for structural dy-

namics with improved numerical dissipation: the generalized-αmethod. Trans.

ASME J. Appl. Mech., 60(2):371–375, 1993.

20 CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS

[Cia88] P. G. Ciarlet. Mathematical Elasticity, Vol. 1. North-Holland, Amsterdam,

1988.

[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.

[CP03] R. W. Clough and J. Penzien. Dynamics of Structures. McGraw-Hill, third

edition, 2003.

[Emm99] E. Emmrich. Discrete versions of Gronwall’s lemma and their application to the

numerical analysis of parabolic problems. Preprint 637, Technische Universit¨at

Berlin, Germany, 1999.

[Eˇ

ST13] E. Emmrich, D. ˇ

Siˇska, and M. Thalhammer. On a full discretisation for non-

linear second-order evolution equations with monotone damping: construction,

convergence, and error estimates. Technical report, University of Liverpool,

2013.

[ET10a] E. Emmrich and M. Thalhammer. Convergence of a time discretisation for

doubly nonlinear evolution equations of second order. Found. Comput. Math.,

10(2):171–190, 2010.

[ET10b] E. Emmrich and M. Thalhammer. Stiffly accurate Runge-Kutta methods for

nonlinear evolution problems governed by a monotone operator. Math. Comp.,

79(270):785–806, 2010.

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

(AMS), Providence, second edition, 1998.

[GC01] M. G´eradin and A. Cardona. Flexible Multibody Dynamics: A Finite Element

Approach. John Wiley, Chichester, 2001.

[GGZ74] H. Gajewski, K. Gr¨oger, and K. Zacharias. Nichtlineare Operatorgleichungen

und Operatordifferential-Gleichungen. Akademie-Verlag, 1974.

[Hug87] T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite

Element Analysis. Dover Publications, 1987.

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

Numerical Solution. European Mathematical Society (EMS), Z¨urich, 2006.

[LMT13] R. Lamour, R. M¨arz, and C. Tischendorf. Differential-algebraic equations: a

projector based analysis. Springer-Verlag, Heidelberg, 2013.

[LP86] P. L¨otstedt and L. R. Petzold. Numerical solution of nonlinear differential equa-

tions with algebraic constraints. I. Convergence results for backward differenti-

ation formulas. Math. Comp., 46(174):491–516, 1986.

[Meh13] V. Mehrmann. Index concepts for differential-algebraic equations. In T. Chan,

W.J. Cook, E. Hairer, J. Hastad, A. Iserles, H.P. Langtangen, C. Le Bris,

P.L. Lions, C. Lubich, A.J. Majda, J. McLaughlin, R.M. Nieminen, J. Oden,

P. Souganidis, and A. Tveito, editors, Encyclopedia of Applied and Computa-

tional Mathematics. Springer-Verlag, Berlin, 2013.

[New59] N. M. Newmark. A method of computation for structural dynamics. Proceedings

of A.S.C.E., 3, 1959.

[Ria08] R. Riaza. Differential-algebraic systems. World Scientific Publishing Co. Pte.

Ltd., Hackensack, 2008.

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

Birkh¨auser Verlag, Basel, 2005.

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

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

CONVERGENCE OF THE ROTHE METHOD IN ELASTODYNAMICS 21

[SGS06] W. Schiehlen, N. Guse, and R. Seifried. Multibody dynamics in computational

mechanics and engineering applications. Comput. Method. Appl. M., 195(41–

43):5509–5522, 2006.

[Sha97] A. A. Shabana. Flexible multibody dynamics: review of past and recent devel-

opments. Multibody Syst. Dyn., 1(2):189–222, 1997.

[Sim98] B. Simeon. DAEs and PDEs in elastic multibody systems. Numer. Algorithms,

19:235–246, 1998.

[Sim00] B. Simeon. Numerische Simulation Gekoppelter Systeme von Partiellen und

Differential-algebraischen Gleichungen der Mehrk¨orperdynamik. VDI Verlag,

D¨usseldorf, 2000.

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

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

[Sim13] B. Simeon. Computational flexible multibody dynamics. A differential-algebraic

approach. Differential-Algebraic Equations Forum. Springer-Verlag, Berlin,

2013.

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

Holland, Amsterdam, 1977.

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

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

Habilitationsschrift, Humboldt-Universit¨at zu Berlin, 2003.

[Tr¨o09] F. Tr¨oltzsch. Optimale Steuerung partieller Differentialgleichungen: Theorie,

Verfahren und Anwendungen. Vieweg+Teubner Verlag, Wiesbaden, 2009.

[Wil98] E. L. Wilson. Three Dimensional Static and Dynamic Analysis of Structures: A

Physical Approach with Emphasis on Earthquake Engineering. Computers and

Structures Inc., Berkeley, 1998.

[Wlo87] J. Wloka. Partial Differential Equations. Cambridge University Press, Cam-

bridge, 1987.

[Zei90] E. Zeidler. Nonlinear Functional Analysis and its Applications IIa: Linear

Monotone Operators. Springer-Verlag, New York, 1990.

∗Institut f¨

ur Mathematik MA4-5, Technische Universit¨

at Berlin, Straße des 17. Juni

136, 10623 Berlin, Germany

E-mail address:[email protected]