Analyzing the stability behaviour of DAE solutions and their approximations

(1)

Analyzing the stability behaviour of DAE solutions and their approximations

Roswitha Marz

Antonio R. Rodrguez Santiesteban and Humboldt-University Berlin

Institute of Mathematics Unter den Linden 6 D-10099 Berlin, Germany

1

(2)

1 Introduction

Usually, lower-index dierential-algebraic equations (DAEs)

f(x⁰(t) x(t) t) = 0 (1.1)

are integrated by numerical algorithms that have been developed originally for regular ordinary dierential equations (ODEs). There are numerous papers justifying this by convergence proofs, order results etc. (cf. 1], 2], 3]). In particular, the backward dierentiation formula (BDF)

f

0

@1 h

k

X

j⁼⁰jxn^;j xn tn

1

A= 0 (1.2)

but also implicit Runge-Kutta methods (IRK) xn=xn^;1+h^X^s

j⁼¹jX_nj⁰ (1.3 a)

f

0

@X_ni⁰ xn^;1+h^X^s

j⁼¹ijX_nj⁰ tn^;1+cih

1

A= 0 i= 1 ::: s (1.3 b) are used in applications with really great success. Step by step each method provides numerical approximations xn of the true solution values x(tn) tn=tn^;1+h. The convergence results mentioned above concern the behaviour of the global error on a compact interval, say t⁰ T], as the stepsize of the discretizationt⁰ < t¹ < <

tn=T tends to zero.

However, what do we know about the error x(tn)^;xn if the stepsize h is constant and n ^! ¹, that is, tn ^! ¹? When integrating regular ODEs numerically, one tries to carefully match the dynamics of the numerical algorithm with the dynamical behaviour of the true solution. How could this be done in case of DAEs?

Classical linear stability theory is concerned with the analysis of approximating the so-called scalar linear test equation

z⁰(t) =z(t) Re <0: (1.4) Basic notions like the region of absolute stability, error growth function etc., rely on (1.4) and apply, via similarity transforms, to constant coecient regular ODEs

x⁰(t)^;Wx(t) = 0: (1.5)

Considering DAEs, the respective constant coecient system

Ax⁰(t) +Bx(t) = 0 (1.6)

3

(4)

has a singular leading coecient matrixA, but a regular matrix pencil^fA B^g, i. e., det(A+B)⁶0. There are nonsingular matricesE F that transform (1.6) into its so-called Kronecker normal form

I J

!

x~⁰(t) +

;W I

!

x~(t) = 0 (1.7)

where EAF =

I J

!

EBF =

;W I

!

x~ = F^;1x and J is a nilpotent block, ind(J) = ind^fA B^g.

Obviously, the solution of the homogeneous equation (1.6) consists of its "dynamical part" only, i. e.,

x(t) =Fx~(t) = F

u~(t) 0

!

where ~u(t) solves the regular ODE ~u⁰(t) =Wu~(t).

On the other hand, we may apply the same transformations E F to decouple the BDF applied to (1.6), that is,

A¹h

k

X

j⁼⁰jxn^;j+Bxn = 0 to

I J

!1 h

k

X

j⁼⁰jx~n^;j+

;W I

!

x~n = 0:

Starting with consistent values x⁰ xk^;1 (such that F^;1xj =

u~j

0

!

j = 0 k) we obtain xn=Fx~n =F

u~n

0

!

where ~un is given by the BDF applied to the regular ODE ~u⁰(t) =Wu~(t), i. e., h1

k

X

j⁼⁰ju~n^;j =Wu~n:

In the same way we may proceed with Runge-Kutta methods. Obviously, the rich world of classical stability theory applies to constant coecient DAEs (1.6) via transformation into Kronecker normal form plus similarity transform.

However, unfortunately, the homogeneous coecient DAE (1.6) does not play a similar role in the DAE analysis as equation (1.5) does in the regular ODE-theory.

DAEs represent a much more complex class of problems. The constant coecient 4

(5)

case (1.6) is only a very poor model, which does not reect important geometric features of DAEs at all.

Positive results concerning long term integrations of index-1 DAEs are reported in 4]. There, the leading nullspace is supposed to remain invariant. Moreover, in 4]

it is stressed that varying subspaces may have its eect.

In 5], the following small example was discussed to show the bad eect of a rotating nullspace. The linear index-1 DAE

^;1 t

0 0

!

x⁰(t) +

^;1 t

^;1 t^;1

!

x(t) = 0 (1.8)

with real parameters ⁶= 0 ⁶= 1 has the solutions x¹(t) = (^;1)^;1(1^;t)x²(t)

x²(t) = exp((^;)t)x²(0):

The nullspace of the leading coecient matrix is N(t) :=ⁿz ²IR² : (^;1)z¹+tz² = 0^o:

For ⁶= 0, this nullspace varies with t. Applying the backward Euler method to (1.8) yields

xn¹ = (^;1)^;1(1^;tn)xn²

xn² = ¹⁺¹⁺_h^hxn^;1²:

Obviously, it holds that xn ^! 0 (n ^! ¹) if and only if ^j1 +h^j < ^j1 +h^j. For h = ^;1 the method is not feasible at all. Further, there is a large region on the (h h)-plane where the true solution x(tn) vanishes asymptotically, but the numerical approximation grows unboundedly.

Note that the spectrum of the pencil ^fA(t) B(t)^gis time-invariant, namely (A(t) B(t)) := ^f²IC : det(A(t) +B(t)) = 0^g=^f;^g:

If we put = 0, we obtain again a constant coecient DAE and everything is ne as expected before.

One could think that the strange asymptotic behaviour is due to the rotations of the nullspace of the leading coecient A(t) or the leading partial Jacobian f_x⁰⁰(x⁰ x t) in (1.1). Since this matrix has a constant nullspace or is constant itself in most applications, we can perhaps consider equation (1.8) to be just an academic example.

However, similar phenomena of a "wrong numerical dynamical behaviour" may arise even in the case of a constant leading coecient matrix. Such phenomena are discussed in 6], where linear index-2 DAEs are considered.

5

(6)

Roughly speaking, in 6] it is shown that the numerical approximations t the dynamical behaviour of the true solution well, supposed both the characteristic subspaces N¹(t) andS¹(t) do not vary with t in fact.

The aim of the present paper is to show that the same invariance condition forS¹(t) will do. This weaker sucient condition ensures also an appropriate reection of the exponential decay. Since e.g. in circuit simulation S¹(t) is often invariant but N¹(t) is not, our new approach is of particular interest for those applications.

This paper is organized as follows. In Section 2 the related analytical background is given. Useful analytical decoupling and reduction techniques are discussed. As a by product, Theorem 2.2 provides new solvability results with respect to weaker smoothness demands. Also contractivity is now discussed under those weaker smoothness conditions.

In Section 3 the analytical decoupling and reduction techniques are applied to BDF and IRK methods, to prove that the time invariance of the subspace S¹(t) will do in fact.

Finally we demonstrate by two characteristic examples that are even in Hessenberg form, how a rotating subspace S¹(t) eects the dynamic reection.

2 Analysis of linear index-2 DAE's

2.1 Fundamentals

Consider the linear equation

A(t)x⁰(t) +B(t)x(t) =q(t) t²J := t⁰ ¹) (2.1) where the coecients are continuous and the matrix A(t) ² L(IR^m) is singular for all t²J, but has constant rankr. Introduce the basic subspaces

N(t) := kerA(t)

S(t) := ^fz ²IR^m :B(t)z ² imA(t)^g:

Obviously, each solution of the homogeneous equation satises x(t)²S(t) t²J:

In the following, we assumeN(t) to be spanned bym^;r continuously dierentiable base functions. Then, there is a C¹ matrix functionQ : J ^!L(IR^m) that projects IR^m pointwise onto N(t) i. e., Q(t)² =Q(t) im Q(t) =N(t) t²J.

In what follows, Q denotes such a C¹ projector function onto N, and P :=I^;Q. Recall (e.g. 2]) the equation

A(t)^f(Px)⁰(t)^;P⁰(t)x(t)^g+B(t)x(t) =q(t) t ²J (2.2) 6

(7)

to be the precise formulation of (2.1). We are looking for continuous functions x( ) :J ^!IR^m that have continuously dierentiable parts (Px)( ), and that satisfy (2.2) pointwisely. Denote this function space by

C_N¹ :=ⁿy²C(J IR^m) :Py²C¹(J IR^m)^o:

Next, we introduce further subspaces which are relevant for index-2 DAEs 2], namely

N¹(t) := kerA¹(t)

S¹(t) := ^fz²IR^m :B(t)P(t)z ² imA¹(t)^g A¹(t) := A(t) + (B(t)^;A(t)P⁰(t))Q(t) t²J:

Denition:

The DAE (2.1) is said to be index-2 tractable if dim(N(t)^\S(t)) = >0 N¹(t)^\S¹(t) =^f0^g t²J.

Recall that index-2 tractability generalizes the case of Kronecker index 2. Let Q¹(t) denote the projector onto N¹(t) alongS¹(t) t²J.

Due to 4], Lemma A.13, the matrix G²(t) :=A¹(t) +B(t)P(t)Q¹(t)

remains nonsingular now and it holds that Q¹(t) = Q¹(t)G²(t)^;1B(t)P(t) on J. Since we may represent

kerA¹(t) = (I^;P(t)A(t)⁺(B(t)^;A(t)P⁰(t))Q(t))(N(t)^\S(t)) we know N¹(t) to have the same constant dimension as N(t)^\S(t).

In the following, we drop the argument t if possible.

By straightforward calculations we prove the relations

G^;1² A =P¹P G^;1² B =G^;1² BPP¹+Q¹+Q+P¹PP⁰Q:

Hence, scaling equation (2.1) by G^;1² leads to

P¹P((Px)⁰^;P⁰x) +G^;1² BPP¹x+Q¹x+Qx+P¹PP⁰Qx=G^;1² q: (2.3) Multiplying (2.3) by PP¹ QP¹ and Q¹, respectively, and carrying out simple com- putations we obtain the decoupled system

PP¹(Px)⁰+PP¹G^;1² BPP¹x = PP¹G^;1² q (2.4)

;QQ¹(Px)⁰+QP¹G^;1² BPP¹x+Qx = QP¹G^;1² q (2.5)

Q¹x = Q¹G^;1² q: (2.6)

7

(8)

The System (2.4), (2.5), (2.6) provides the basic idea of what the solutions of the DAE (2.1) look like in case of PP¹ PQ¹ and PQ¹G^;1² q beingC¹. Then (2.4), (2.5) may be rewritten as

(PP¹x)⁰ ^;(PP¹)⁰PP¹x+PP¹G^;1² BPP¹x=PP¹G^;1² q+ (PP¹)⁰PQ¹G^;1² q (2.4') Qx=^;(QP¹G^;1² B+QQ¹(PQ¹)⁰)PP¹x + QQ¹(PQ¹G^;1² q)⁰+QP¹G^;1² q

; QQ¹(PQ¹)⁰PQ¹G^;1² q (2.5') and the solution is composed by x = PP¹x+PQ¹x+Qx where the component PP¹xsolves the inherent regular ODE (2.4'),Qx and PQ¹x are given by (2.5') and (2.6), respectively. (cf. 2]).

In the next section, we will generalize the solvability results given in 8] in the sense that we will do with lower smoothness. The new inherent regular ODE, which uses a dierent projector instead of PP¹, permits to obtain new insights into the asymptotic stability behaviour of numerical approximations (Section 3).

Note that, by construction, P(t)P¹(t) projects onto P(t)S¹(t) along N(t)N¹(t), i. e., it is related to the decomposition

IR^m =P(t)S¹(t)N(t)N¹(t): (2.7)

2.2 Solvability and decoupling

Now we also use the orthoprojectorV(t), which projects alongS¹(t). Then,I^;V(t) is the orthoprojector ontoS¹(t). Because ofN(t) S¹(t), it holds thatV(t)Q(t) = 0, and

(t) :=P(t)(I^;V(t))

is a projector, too. Obviously, we have im (t) = P(t)S¹(t). Since the two projectors (t) and P(t)P¹(t) have the same image space P(t)S¹(t), it holds that PP¹ = PP¹ =PP¹.

Additionally, we denote the orthoprojector onto im A¹(t) by I^;W(t).

As we will realize below, the projection W extricates exactly that derivative free part of equations that has to be dierentiated once when reducing the index.

The following example of a Hessenberg form DAE is to make the meaning of the dierent projectors more transparent.

Example:

Given a Hessenberg form DAE of size 2 that contains n¹ equations with derivatives and n² derivative-free ones

x⁰¹+ B¹¹x¹+B¹²x² = q¹ B²¹x¹ = q²

)

B²¹B¹² nonsingular.

8

(9)

Here we have a constant nullspace N N =^fz :z¹ = 0^g further Q=

0 0 0 I

!

A¹ =

I B¹² 0 0

!

W =

0 0 0 I

!

S(t) =S¹(t) =^fz :B²¹(t)z¹ = 0^g S(t)^\N =N N¹(t) =^fz:z¹ +B¹²(t)z² = 0^g

N¹(t)^\S¹(t) =^fz :z¹+B¹²(t)z² = 0 B²¹(t)z¹ = 0^g=^f0^g due to the nonsingularity of B²¹B¹².

The projector PP¹ is given by PP¹ =

I^;L 0

0 0

!

where L =B¹²(B²¹B¹²)^;1B²¹ is also a projector (L(t) projects IRⁿ¹ onto im B¹²(t) along kerB²¹(t)). The orthoprojector along S¹(t) is now

V(t) =

B²¹(t)⁺B²¹(t) 0

0 0

!

and (t) =

I^;B²¹(t)⁺B²¹(t) 0

0 0

!

im (t) = P(t)S¹(t) =^fz :B²¹(t)z¹ = 0 z² = 0^g

= kerB²¹(t)^f0^g:

It is well-known that, in Hessenberg systems, all derivative-free equation should be dierentiated when deriving a solution via a reduction step.

Consequently, it seems to be more natural to assume B²¹ ² C¹ (hence, V ²C¹) instead of L²C¹ (i. e., PP¹ ²C¹).

2

Next, we collect some nice properties of our projectors and matrices to be used below:

1) im A¹ = kerW im A¹ = kerG²PQ¹G^;1² , hence there are two projectors W and G²PQ¹G^;1² along imA¹. Therefore

W =WG²PQ¹G^;1² G²PQ¹G^;1² =G²PQ¹G^;1² W:

9

(10)

2) FromWB =WG²PQ¹G^;1² B =WG²PQ¹G^;1² BPQ¹ =WBPQ¹ and PQ¹G^;1² WB =PQ¹ it follows that

kerWB = kerPQ¹ = kerQ¹ =S¹ imWB = im W:

3) V = (WB)⁺WB W = (I^;A¹A⁺¹) = (PQ¹G^;1² )⁺PQ¹G^;1² WB = (PQ¹G^;1² )⁺PQ¹:

4) V =V PQ¹ PQ¹ =PQ¹V V =V P V =V Q¹ V G^;1² = (WB)⁺WBG^;1² = (WB)⁺WBPQ¹G^;1²

= (WB⁺WG²PQ¹G^;1² BPQ¹G^;1² = (WB)⁺WG²PQ¹G^;1² = (WB)⁺W. The third property indicates that, for continuously dierentiablePQ¹G^;1² andPQ¹, we also have continuously dierentiable W and V. However, as we could realize by means of the special Hessenberg form DAE above, the opposite is not true. In this sense, Theorem 2.2 below generalizes the related results of 8].

To get a better insight we reconsider the decoupled form (2.4), (2.5), (2.6) of the DAE (2.1) and try to reformulate it in such a way that (2.4) returns into a regular ODE for the solution component x. Now, we do not assume PP¹ PQ¹ to be from C¹ as we did in order to obtain (2.4'), (2.5'), but we allow only P and V to be continuously dierentiable now.

Equation (2.6) yields immediately

V x= (WB)⁺Wq: (2.8)

To realize what equation (2.4) looks like in more detail we derive PP¹(Px)⁰ = PP¹( +PV)(Px)⁰ = (Px)⁰+PP¹V(Px)⁰

= (x)⁰^;⁰Px+PP¹(V x)⁰^;PP¹V⁰Px:

Consequently, (2.4) can be rewritten as

(x)⁰^;⁰(x+ PV x)^;PP¹V⁰(x+PV x) +PP¹(V x)⁰

+ PP¹G^;1² BPP¹(x+PV x) =PP¹G^;1² q: (2.9) Analogously, with QQ¹(Px)⁰ = QQ¹V(Px)⁰ = QQ¹(V x)⁰ ^; QQ¹V⁰(x+ PV x), equation (2.5) can be reformulated:

;QQ¹(V x)⁰+QQ¹V⁰(x+PV x)+QP¹G^;1² BPP¹(x+PV x)+Qx=QP¹G^;1² q:(2.10)

10

(11)

If x²C_N¹ solves the DAE (2.1), then (WB)⁺Wq=V x=V Pxbelongs toC¹ since V and Px do so. We are allowed to replace V x and (V x)⁰ in (2.9) and (2.10) by means of (2.8). This gives rise to consider the ODE

u⁰^;⁰u ^;PP¹V⁰u+PP¹G^;1² BPP¹u

= ⁰P(WB)⁺Wq+PP¹V⁰P(WB)⁺Wq^;PP¹G^;1² BPP¹(WB)⁺Wq

;PP¹((WB)⁺Wq)⁰+PP¹G^;1² q (2.11) and to call it an inherent regular ODE. If we multiply equation (2.11) by (I^;), we obtain

(I ^;)u⁰ ^;(I ^;)⁰u= 0 hence

((I^;)u)⁰+ ⁰(I^;)u= 0:

It results that if u ²C¹ solves (2.11), then the function ! = (I^;)u satises the homogeneous regular ODE !⁰+⁰!= 0. If, additionally,!(t) = (I^;(t))u(t) = 0, then !(t) vanishes identically, i. e., u(t) = (t)u(t) t ² J. This proves the following assertion.

Lemma 2.1

The subspace im (t) IR^m represents an invariant subspace of the regular inherent ODE (2.11).

As far as the homogeneous case q = 0 is concerned, we know from (2.8) that each solution of Ax⁰ +Bx = 0 has a trivial component V x = 0. Using (2.9), (2.10) and taking into account the inherent regular ODE we may express each solution as x = x +Qx = canu, where u solves (2.11) as well as the initial condition u(to)² im (t⁰), and

can:=KPP¹ K := (I^;QP¹G^;1² BPP¹^;QQ¹V⁰PP¹):

The matrix function K is nonsingular, thus ker can(t) = kerP(t)P¹(t) t ²J. It is easily checked that can is also a projector. canis said to be the canonical projector for the index-2 case. It projects onto the solution space of the homogeneous DAE along N N¹, as we will see below.

Theorem 2.2

Given an index-2 DAE (2.1) with continuously dierentiable subspaces N S¹ and im A¹. Additionally, let WB be C¹ and q²^fp²C:Wp²C¹^g.

(i) Then, the IVP for (2.1) with the initial condition (t⁰)(x(t⁰)^;x⁰) = 0 x⁰ ²IR^m, has exactly one C_N¹-solution.

(ii) Exactly one solution of the homogeneous equation at t⁰ passes through each x⁰ ²im can(t⁰) .

11

(12)

Proof:

(i) Denote by u ² C¹ the solution of the regular ODE (2.12) that satises the initial conditionu(t⁰) = (t⁰)x⁰. Then we compose the functionx=u+v+w v := (WB)⁺Wq²C¹

w:=QP¹G^;1² q+QQ¹v⁰^;QQ¹V⁰(u+Pv)^;QP¹G^;1² BPP¹(u+Pv): Due to this construction we have v =V v w=Qw Px=u+Pv ²C¹ x²C (t⁰)x(t⁰) =u(t⁰) = (t⁰)x⁰.

Finally, straightforward checking shows thatxsatises the DAE (2.1), indeed.

(ii) Solve the special IVP for q = 0 and x⁰ := x⁰ = can(t⁰)x⁰. Its solution is x= canu, and we have, in particular,

x(t⁰) = can(t⁰)u(t⁰) = can(t⁰)(t⁰)x⁰

= can(t⁰)(t⁰)can(t⁰)x⁰ = can(t⁰)(t⁰)(PP¹)(t⁰)x⁰

= can(t⁰)(PP¹)(t⁰)x⁰

= can(t⁰)x⁰ =x⁰:

2

Corollary:

On each compact interval t⁰ T], the perturbation index of an index-2 tractable DAE is two.

Proof:

For each T > t⁰, there is a constantKT such that the estimation

kx^kT KT

n

j(t⁰)x⁰^j+^kq^kT +^k(V q)⁰^kT

o (2.12)

holds for all IVP solutions with x⁰ ²IR^m q²C V q²C¹, where

ky^kT := max^fjy(t)^j:t²t⁰ T]^g.

2

Example

: For the Hessenberg form DAE x⁰¹+ B¹¹x¹+B¹²x² = q¹

B²¹x¹ = q²

)

the conditions of Theorem(2.2) are satised if q¹ ²C q² ²C¹ and B²¹²C¹. The initial condition reads in detail

0 = (t⁰)(x(t⁰)^;x⁰) =

(I^;B²¹(t⁰)⁺B²¹(t⁰))(x¹(t⁰)^;x⁰¹) 0

!

:

2.3 Index reduction by dierentiation

Considering the DAE (2.1) we observe immediately that each solution has always to proceed in the set

M(t) := ^fz ²IR^m :B(t)z^;q(t)² imA(t)^g=

fz ²IR^m : (I^;W(t))(B(t)z^;q(t))² imA(t) W(t)(B(t)z^;q(t)) = 0^g 12

(13)

which is called a constraint manifold. Under the conditions of Theorem 2.2, the part of equations described by

WBx=Wq

can be dierentiated to obtain

WBx⁰+ (WB)⁰x= (Wq)⁰:

Here, WBx⁰ stands forWB^f(Px)⁰^;P⁰x^g:On the other hand, due to (2.1) we may express Px⁰ := (Px)⁰ ^;P⁰x=PA⁺(q^;Bx): In consequence, the relation

WBA⁺(q^;Bx) +W(WB)⁰x=W(Wq)⁰

has to be satised additionally, i. e., the solution has also to proceed in the so-called hidden constraint manifold

H(t) := ^fz ²IR^m : (^;W(t)B(t)A(t)⁺B(t) +W(WB)⁰(t))z

=W(Wq)⁰(t)^;W(t)B(t)A(t)⁺q(t)^g: As we will see below, the set

M

1(t) := ^M(t)^\^H(t)

is characteristic of the DAE (2.1). ^M¹ contains all solutions (for given q), but it is lled by those solutions. One could call^M¹(t) the state manifold. For the case of a homogeneous equation, S(t) = ^M(t) is trivially given. Moreover, there is a closed relationship between our subspace S¹(t) and the state manifold^M¹(t).

Lemma 2.3

For a homogeneous index-2 tractable DAE (2.1) with continuously dif- ferentiable N S¹ im A¹ WB it holds that

M

1(t)^jq⁼⁰ =ⁿz ²S¹(t) :Q(t)z =^;((QP¹G^;1² B+QQ¹G^;1² (WB)⁰)PP¹)(t)z^o and

P(t)^M¹(t)^jq⁼⁰ =P(t)S¹(t): (2.13)

Proof:

The second relation is a simple consequence of the representation of ^M¹(t), which we want to realize now. We drop the argument t.

z ²^M¹ means by denition:

;WBA⁺Bz+W(WB)⁰z = 0 Bz+Aw= 0 for a certain w=Pw. This is equivalent to

WBw+W(WB)⁰z = 0 w=^;PA⁺Bz PQ¹z = 0 Qz^;QQ¹w+QP¹G^;1² BPP¹z+QP¹PP⁰Qz = 0, and to

13

(14)

PQ¹z = 0 w=^;PA⁺Bz 0 =WBw+W(WB)⁰PP¹z+W(WB)⁰Qz = WBw+W(WB)⁰PP¹z+WBP⁰Qz

Qz^;QQ¹w^;QQ¹P⁰Qz+QP¹G^;1² BPP¹z = 0:

Because of QQ¹G^;1² WB = QQ¹G^;1² B = QQ¹ QQ¹G^;1² W = QQ¹G^;1² , we nd z ²^M¹ to be characterized by

PQ¹z = 0 w =^;PA⁺Bz

QQ¹w+QQ¹P⁰Qz+QQ¹G^;1² (WB)⁰PP¹z= 0 Qz = QQ¹w+QQ¹P⁰Qz^;QP¹G^;1² BPP¹z

= ^;QQ¹P⁰Qz+QQ¹P⁰Qz^;QQ¹G^;1² (WB)⁰PP¹z^;QP¹G^;1² BPP¹z

= ^;QQ¹G^;1² (WB)⁰PP¹z^;QP¹G^;1² BPP¹z

2

Note that the linear space ^M¹(t)^jq⁼⁰ represents the tangent space for ^M¹(t). Of course, this is much more important for nonlinear problems. Lemma 2.3 shows P(t)S¹(t) to be a kind of practical substitute of the tangent space. It should be stressed once more that we are looking forC_N¹ solutions and that there is no natural need for Q-components of the elements of the tangent spaces.

Next, rewrite the DAE (2.1) as

Ax⁰+ (I^;W)(Bx^;q) +W(Bx^;q) = 0

and replace the part W(Bx^;q) = 0 by its dierentiated form W^fWBx⁰+ (WB)⁰x^;(Wq)⁰^g= 0

to compose the new DAE

(A+WB)x⁰+ ((I^;W)B+W(WB)⁰)x= (I ^;W)q+W(Wq)⁰: (2.14) Denote ~A := A +WB B^~ := (I ^; W)B +W(WB)⁰ and consider (2.14) in more detail. Since (A(t) +W(t)B(t))z = 0 decomposes into Az = 0 WBz = 0, we have ker ~A(t) kerA(t) = N(t). In consequence, the solutions of (2.15) belong to the same class C_N¹ as the solutions of (2.1). Further, we have

S~(t) =^fz ²IR^m: ~B(t)z ² im ~A(t)^g

=^fz ²IR^m: (I^;W(t))B(t)z ² imA(t)

W(t)(WB)⁰(t)z =W(t)B(t)A(t)⁺(I^;W(t))B(t)z^g

f

M(t) =^fz ²IR^m : (I ^;W(t)(B(t)z^;q(t))² im A(t)

W(t)(WB)⁰(t)z^;W(t)(Wq)⁰(t) = W(t)B(t)A(t)⁺(I^;W(t))(B(t)z^;q(t))^g

f

M(t)^\^fz²IR^m :W(t)(B(t)z^;q(t)) = 0^g=^M¹(t): 14

(15)

Theorem 2.4

Given the conditions of Theorem 2.2.

(i) Then equation (2.14) is index-1 tractable.

(ii) Exactly one solution at t ² J passes through each x ²^M^f(t) . (iii) ^M¹(t) and ^M(t) form invariant subsets for (2.14).

Proof:

(i) We check the nonsingularity of the matrix ~G¹(t) ²L(IR^m) ~G¹ :=A+WB + ((I^;W)B+W(WB)⁰)Q:Again, we drop the t-argument. ~G¹z = 0 yields Az+(I^;W)BQz = 0 WBz =^;W(WB)⁰Qz =^;WBP⁰Qz, hencePQ¹z =

;PQ¹P⁰Qz (A+BQ)z = 0.

Because of (A+BQ)(I^;PP⁰Q) =A¹ A¹(I+PP⁰Q) =A+BQ, the element (I +PP⁰Q)z = ~z belongs to kerA¹ =N¹ , i. e., ~z =Q¹z:~ On the other hand, PQ¹z +PQ¹P⁰Qz = 0 implies PQ¹(I +PP⁰Q)z = 0, i. e., PQ¹z~ = 0 and nally ~z = 0 z = 0:

(ii) Due to the index-1 property, this assertion is now given by the index-1 results, e. g. in 2].

(iii) Let x denote a solution of the index-1 tractable DAE (2.14) and x(t)²^M(t): Then, (2.14) gives

Ax⁰+ (I^;W)Bx= (I^;W)q WBx⁰+W(WB)⁰x=W(Wq)⁰: Consider the function :=W(Bx^;q). Derive

⁰ = (WB)⁰x+WBx⁰ ^;(Wq)⁰

= (WB)⁰x^;W(WB)⁰x+W(Wq)⁰^;(Wq)⁰

= (I^;W)(WB)⁰x^;(I^;W)(Wq)⁰

= ^;(I^;W)⁰(WBx^;Wq) =W⁰:

Since x(t)² ^M(t) implies(t) =W(t)(B(t)x(t)^;q(t)) = 0, the function vanishes identically, hence Ax⁰+Bx=q is satised.

Therefore, x(t) ² ^M(t) implies x(t) ² ^M(t) but also x(t) ² ^M¹(t) for all t ²J.

2

15

(16)

Corollary:

The index-2 tractable DAE has dierentiation index two.

Example:

For the Hessenberg system x⁰¹ + B¹¹x¹ + B¹²x² = q¹

B²¹x¹ = q²

)

we have now

M(t) = ^fz :B²¹(t)z¹ ^;q²(t) = 0^g

H(t) = ^fz :^;B²¹(t)B¹¹(t)z¹^;B²¹(t)B¹²(t)z²+ +B²¹⁰ (t)z¹ =q⁰²(t)^;B²¹(t)q¹(t)^g

M

1(t) = ^fz :B²¹(t)z¹ =q²(t) (B²¹(t)B¹²(t))z² =

B²¹⁰ (t)z¹ ^;B²¹(t)B¹¹(t)z¹^;q²⁰(t) +B²¹(t)q¹(t)^g

~

M(t) = ^fz :B²¹⁰ (t)z¹ ^;q²⁰(t) =B²¹(t)(B¹¹(t)z¹ +B¹²(t)z²^;q¹(t))^g

~

M(t) ^\ ^fz :B²¹(t)z¹ =q²(t)^g=^M¹(t): The index-reduced equation (2.15) is, as expected,

x⁰¹ + B¹¹x¹ + B¹²x² = q¹ B²¹x⁰¹ + B²¹⁰ x¹ = q⁰²

)

2.4 Contractivity

Now we turn to the asymptotic behaviour of the solutions of the homogeneous equations

Ax⁰+Bx= 0: (2.15)

As we know, the solutions may be represented by x=KPP¹u= canu

where u solves the inherent regular ODE

u⁰ = ⁰u+PP¹V⁰u^;PP¹G^;1² BPP¹u (2.16) u(t⁰)² im (t⁰):

Recall once more that im (t) is an invariant subspace of the ODE (2.16), i. e., u(t) = (t)u(t) t²J:

16

(17)

Obviously, the stability behaviour of x is mainly governed by the dynamics of the inherent regular ODE, that is, by u. Hence, it would be nice to formulate criteria on how the ow of regular ODE behaves within an given invariant subspace.

The matrix coecient of the inherent regular ODE (2.16) has to be expected to be- come singular. In particular, for constant subspacesN andS¹, it is simply the matrix PP¹G^;1² BPP¹, which has at least the nullspaceNN¹ of dimensionm^;r+ 2.

The standard approaches to estimate how the solutions grow usually relate to all solutions of (2.16), included those solutions that do not start in im (t⁰). No exponential decay can be realized in this way. The standard techniques are somewhat coarse. Hence, we try to improve them by relating all things to the invariant subspaces.

For a given subspace U IR^m we introduce the matrix-semi-norm

kG^k^U := max

(

jGz^j

jz^j ^:z ⁶= 0 z ²U

)

such that ^jGz^j^kG^k^U^jz^j holds for allz ²U.

Then, we dene the logarithmic matrix norm relative to the subspace U by ^U(G) := lim_h

!0

h1⁽^kI+hG^k^U ^;1): (2.17)

This logarithmic norm relative toU is well-dened and has similar properties as the standard version of Dahlquist for U = IR^m (cf.(9]).In fact this generalization is in fact very natural and straightforward. On the other hand, it is worth mentioning that ^U(G) appears to be a special case of the logarithmic norm proposed in 10]

for matrix pencils.

Lemma 2.5

Given a regular linear ODE

u⁰(t) +M(t)u(t) = 0 t ²J = t⁰ ¹)

which has the invariant subspace U(t) IR^m t²J. Let the function (t) :=^Z^t

t⁰ ^U⁽^s⁾(^;M(s))ds t²J be well-dened, i. e., the integrals do exist.

Then, for all ODE solutions starting with initial values u(t⁰)²U(t⁰), the inequality

ju(t)^je⁽^t⁾^ju(t⁰)^j tt⁰ is valid.

17

(18)

Proof:

This proof follows the lines of the standard theory 9], too. Given a solution with u(t⁰)²U(t⁰) we have, u(t)²U(t) for allt²t⁰ ¹).

Denote m(t) := ^ju(t)^j t ²t⁰ ¹), and derive

m(t+h) = ^ju(t) +hu⁰(t) +o(h)^j=^j(I^;hM(t))u(t) +o(h)^j

k(I^;hM(t)^k^U⁽^t⁾^ju(t)^j+o(h) thus

h1⁽m(t+h)^;m(t)) 1

h⁽^kI^;hM(t)^k^U⁽^t⁾^;1)m(t)^;o(h⁰): Letting h^!0, we obtain

D⁺m(t)^U⁽^t⁾(^;M(t))m(t) t ²t⁰ ¹)

where D⁺m(t) denotes the respective Dini derivative. By Peano's Lemma it follows immediately that

m(t)e⁽^t⁾m(t⁰):

2

Theorem 2.6

Let (2.15) be index-2 tractable with N S¹ and WB from C¹. Then, the estimation

jx(t)^j^kcan(t)^ke⁽^t⁾^j(t⁰)x(t⁰)^j tt⁰ (2.18) with (t) := _t^R^t

0

⁽^PS¹⁾⁽⁾(^;M())d t t⁰, holds true for each solution, provided that (t) is well-dened. Here,

M :=^;⁰^;PP¹V⁰+PP¹G^;1² BPP¹

is the coecient matrix of the inherent regular ODE (2.16).

Proof:

Applying Lemma 2.5 to (2.16), which has the invariant subspace im (t) = P(t)S¹(t), yields ^ju(t)^je⁽^t⁾^ju(t⁰)^j tt⁰, provided that u(t⁰)² im (t⁰).

The DAE solutions have the representationx(t) = can(t)u(t), andu(t⁰) = (t⁰)x(t⁰).

2

18

(19)

Corollary:

If there is a > 0 such that ⁽^PS¹⁾⁽^t⁾(^;M(t)) ^; t t⁰, all DAE solutions satisfy the inequality

jx(t)^j^kcan(t)^ke^;⁽^t^;^t⁰⁾^j(t⁰)x(t⁰^j t t⁰: (2.19) Hence, if can(t) grows "moderately" like polynomials or is bounded by

kcan(t)^kCe^t tt⁰, with < , then the solutions decrease exponentially.

To realize whether numerical approximations reect the stability behaviour of the true solution well, one often uses the notion of contractivity and so-called one-sided Lipschitz conditions (cf. (9]). This approach applies also to DAEs provided that we relate the things again to our invariant subspaces.

Standard numerical integration methods applied to index-2 DAE are expected to work well if the basic nullspace N is constant. However, it is well-known that these methods may fail in case this nullspace rotates with time. This is why we suppose P⁰ = 0 now.

If the vector norm used to dene the logarithmic norm is related to an inner product, it holds that

hGz zⁱ^U(G)^jz^j² for allz ²U:

Denition:

An index-2 tractable DAE with constant P is said to be contractive if there are an inner product ^h ⁱ and a constant >0 such that the inequality

hy Pxⁱ^;^jPx^j² (2.20)

is valid for all y x²IR^m t²t⁰ ¹), with

A(t)y+B(t)x= 0 Qy= 0 V(t)y=V(t)⁰(t)x: (2.21) This denition is given in 6] for the case of a continuously dierentiable Q¹ and a real scalar product. If Q¹ ²C¹, we may make use of

Q¹(t)y =Q¹(t)V(t)y =Q¹(t)⁰(t)x=^;Q⁰¹(t)(t)x=^;Q⁰¹(t)PP¹(t)x and arrive at the same expression as used in 6] instead of (2.21).

By decoupling A(t)y+B(t)x= 0, we nd that (2.21) leads to y+M(t)(t)x= 0 V(t)x= 0 Px= (t)x:

Therefore, (2.20) reads then

h;M(t)(t)x (t)xⁱ^;^j(t)x^j²

which means a contractivity condition for the inherent regular ODE relative to the invariant subspace im (t) =PS¹(t).

19

(20)

3 Analyzing numerical integration methods for linear index-2 DAEs

In this section we derive conditions for preserving stability properties of numerical methods when solving index-2 linear DAEs. Let the conditions of Theorem 2.2 be fullled in the following. As in 6] we assume N(t) to be constant, otherwise it is well known that those methods will fail. We shall consider the discrete version of the decoupling of the last section for two of the most important families of methods.

An analogous analysis was made in 6] using the standard approach mentioned in Section 2.1. In that paper they obtained as a sucient condition for preserving the stability properties that both subspaces N¹ and S¹ have to be time invariant. By the discrete decoupling with the new approach described in Section 2.2, we shall get a weaker condition as in 6]. Namely, the time invariance of S¹ will do.

3.1 BDF

Consider the homogeneous system

A(t)x⁰(t) +B(t)x(t) = 0 t²J := t⁰ ¹): (3.1) For homogeneous index-2 tractable DAEs (3.1), the solutions are given by the expression

x= (I^;QP¹G^;1² BPP¹^;QQ¹V⁰PP¹)PP¹u= canu (3.2) where u solves the regular ODE

u⁰ ^;⁰u^;PP¹V⁰u+PP¹G^;1² BPP¹u= 0 (3.3) with u(t⁰)²Im(t⁰). Suppose P⁰ = 0.

If the projector V is time invariant, or equivalently, = P(I^;V) is so, this solution representation simplies to

x= (I^;QP¹G^;1² BPP¹)PP¹u (3.4)

u⁰+PP¹G^;1² BPP¹u= 0: (3.5)

Given the step size h >0,ti =t⁰+ih, i²IN. The BDF applied to (3.1) reads Ai^Xk

j⁼⁰jxi^;j +hBixi = 0 ik (3.6)

where the starting values x⁰ ::: xk^;1 are supposed to be known.

First, we multiply by ((WB)⁺W)i and obtain

Vixi = 0 ik (3.7)

20

Analyzing the stability behaviour of DAE solutions and their approximations