Managing the drift-off in numerical index-2 differential algebraic equations by projected defect corrections

Managing the drift-o in numerical index-2 dierential algebraic equations by projected defect corrections

Roswitha Marz, Berlin

Abstract

When integrating index-2 dierential-algebraic equations, the given constraint may be failed to be met due to the integration method itself and also due to numer- ical defects in the realization. This so-called drift-o gives rise to bad instabilities.

In 1991 Ascher and Petzold proposed to manage the drift-o caused by symmetric implicit Runge-Kutta methods in Hessenberg systems by means of backprojections onto the constraint. In the present paper, this nice idea is generalized and analyzed in some detail for general index-2 dierential-algebraic equations and, in particular, for quasilinear equations ^a(^{x t})^x0+^g(^{x t}) = 0, as they arise in applications. Now the constraint under consideration is only implicitly given and the backprojection turns out to be rather a projected defect correction.

1 Introduction

In order to apply symmetric discretization problems reasonably also in case of index-2 dierential-algebraic equations (DAEs)

x⁰¹(t) +b¹(x¹(t) x²(t) t) = 0 (1.1)

b²(x¹(t) t) = 0 (1.2)

Ascher and Petzold (1991) proposed to complete standard methods by an additional back- ward projection onto the constraint given by equation (1.2). An already known approxi- mationxl¹ xl² of the true solution valuex¹(tl) x²(tl) is replaced by a new approximation x^newl¹ x^newl² such that the condition

b²(x^new_l¹ tl) = 0 (1.3)

is satised. This is accomplished by means of the ansatz x^new_l¹ = xl¹+ @

@x²b¹ x^new_l¹ x^new_l² tl

(1.4)

x^new_l² = xl²: (1.5)

The defect b²(xl¹ tl) represents the drift o from the constraint manifold given by equa- tion (1.2). Clearly, with the backprojection one aims at reducing this drift o.

The index-2 property of (1.1),(1.2) guarantees that and x^new_l are locally uniquely de- termined by (1.3)-(1.5) (cf. Hairer and Wanner (1991), p. 513) provided that the initial

1

approximations are suciently accurate.

One step of the Newton method applied to the nonlinear system (1.3)-(1.5) for x^new_l starting with the initial guess x^new(0)_l :=xl ⁽⁰⁾ := 0 gives

x^new(1)_l¹ =xl^1;Flb²(xl¹ tl) (1.6) where we denote ^@b_@x¹² _@x^@b21_@x^@b12;1(xl¹ xl² tl) =: Fl. Furthermore, if we take into account that Fl @b@x²¹(xl¹ tl) =: Hl is a projector, then (1.6) turns out to be a kind of projected defect correction. Namely, sinceHlFl =Fl holds true, formula (1.6) means in more detail

Hlx^new(1)_l¹ = Hlxl^1;Fl b²(xl¹ tl) (I^;Hl)x^new(1)_l¹ = (I^;Hl)xl¹:

Thus, only the particular component related to the projectorHl is aected.

The resulting projected Runge-Kutta methods have proved their value in lots of cases.

The idea of a backward projection has also been used for the numerical integration of regular ordinary dierential equations (ODEs) and index-1 DAEs with great success for maintaining given invariants numerically (cf. Shampine (1986), Schulz, Bock and Stein- bach (1996), Eich (1992)).

In Gear (1986), it was shown that the overdetermined system of a regular ODE with an invariant

x⁰¹(t) + (x¹(t) t) = 0 (x¹ (t)) = 0

is equivalent, in some sense, to the special index-2 DAE x⁰¹(t) + (x¹(t) t) + @

@x¹⁽x¹(t))^Tx²(t) = 0 (x¹(t)) = 0

which makes the close relationship between the projected integration method for index-2 DAEs and the backprojection onto invariants more transparent.

The present paper pursues two dierent objectives. On the one hand, the procedure (1.6) of the projected defect correction shall be applied to DAEs that are not of the special Hessenberg form (1.1), (1.2). For more general index-2 DAEs, like those occurring in the circuit simulation when using the charge oriented modied nodal analysis (e.g. Marz and Tischendorf (1996)), the structure is not given as explicitly as in case of Hessenberg systems. The component that has to be improved and which is an analogue to Hlxl¹, as well as the critical part of the defect are implicitly contained in the system and have to be found out rst.

For quasilinear index-2 DAEs

a(x(t) t)x⁰(t) +b(x(t) t) = 0 (1.7) for instance, we propose and investigate a generalization of the correction formulas (1.5), (1.6) that looks like

x^new(1)_l :=xl^;PQ¹l G^;12_l b(xl tl) (1.8) 2

with a scaling matrixG²l and a special projector PQ¹l.

The termPQ¹lG^;12_l b(xl tl) turns out to represent a critical particular component of the defect caused by the given approximationxl in the derivative free part of the DAE (1.7).

In other words, it represents the drift o from a certain constraint implicitly contained in (1.7).

Recall once more that the remarkable thing about Hessenberg form DAEs (1.1), (1.2) is just the explicitly given relevant constraint (1.2).

There are well-known integration methods like the BDF that satisfy a priori the constraint (1.2). Thus, theoretically, applying e.g. the BDF to integrate (1.1), (1.2), there is no drift o caused by the integration method. However, in practical computations we are again confronted with nontrivial defects. Now they represent the residuals of the Newton iterations when solving the nonlinear systems arising per integration step.

The accuracy of numerical approximations generated e.g. by a BDF-code is known to be seriously inuenced by the critical defect components amplied by h^;1_l , where hl denotes the current stepsize. This is why the nonlinear equations arising per integration step have to be solved suciently accurately. Whether an integration code applied to an index-2 DAE performs reliably or not depends strongly on the tolerances for the Newton iterations. The code fails in both cases if tolerances are too precise or too coarse.

The second objective in the present paper is to treat the backprojection as a usual defect correction step within the more critical subsystem of the nonlinear equation system to be solved. As we know, this subsystem is given only implicitly, but can be gured out by projections. Solving this subsystem more precisely by means of defect corrections, we may expect the code to perform well also if the tolerances for the whole nonlinear system are weakened. In the consequence, the integration code gains more reliability and robustness.

The paper is organized as follows:

In Section 2 we develop a projected defect correction method for general implicit index-2 DAE's and show in what sense (1.3)-(1.5) and (1.6) are thus generalized. A detailed error analysis is presented in Section 3, where particularly simple relations result for quasilinear DAEs (1.7).

In Section 4 we describe a method for the numerical calculation of the special projector Q¹l.

Section 5 deals with relations between the projected defect correction and the selective projected methods proposed in Ascher and Spiteri (1994) for the case of semi-explicit index-2 systems.

An experiment in applying the projected defect correction for the improvement of the reliability of the BDF is reported in Section 6.

No daubt, the proposed projected defect corrections allow, as in the case of Hessenberg systems, the use of projected Runge-Kutta methods etc. without having to transform the DAE itself anytime. Here, the eect cannot but be the same as already described in Ascher and Petzold (1991).

3

2 A general projected defect correction formula

Consider the DAE

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

wheref :^Gf^If IR^mIR^mIR^!IR^mis continuous together with its partial Jacobians f_x⁰⁰ f_x⁰ :^Gf ^If ^!L(IR^m):

The nullspace of the leading Jacobian kerf_x⁰⁰(x⁰ x t) is assumed to be constant. Denote N := kerf_x⁰⁰(x⁰ x t)

and introduce projector matricesQ P ²L(IR^m) such thatQ² =Q,imQ =N,P :=I^;Q. In case of the Hessenberg system (1.1), (1.2) the leading nullspace N is simply

N =z¹ z²

2IR^m :z¹ = 0 further we may chooseQ=diag(0 I).

Letx(:) :^{I !}IR^m denote the solution of (2.1) which is to be approximated. The DAE (2.1) is supposed to be index-2 tractable in a certain neighbourhood^GI of the graph

; :=^f((Px)⁰(t) x(t) t) :t^2Ig:

More precisely, this means by denition (e.g. Marz (1995)) that the intersection N ^\ S(x⁰ x t) has a constant but nonzero dimension on ^{G I}, and the subspaces N¹(x⁰ x t) S¹(x⁰ x t) intersect trivially on ^GI, i.e.,

N¹(x⁰ x t)^\S¹(x⁰ x t) =^f0^g (x⁰ x t)^2GI: (2.2) There, the following denotations are used:

S(x⁰ x t) := ^fz ²IR^m :f_x⁰(x⁰ x t)z ²im f_x⁰⁰(x⁰ x t)^g S¹(x⁰ x t) := ^fz ²IR^m :f_x⁰(x⁰ x t)Pz ²im G¹(x⁰ x t)^g G¹(x⁰ x t) := f_x⁰⁰(x⁰ x t) +f_x⁰(x⁰ x t)Q

N¹(x⁰ x t) := kerG¹(x⁰ x t):

Recall that the nullspace N¹(x⁰ x t) has the same dimension as the intersection N ^\ S(x⁰ x t). Moreover, there is a uniquely determined projector matrix function Q¹ : ^G

I !L(IR^m) such that Q¹(x⁰ x t) projects onto N¹(x⁰ x t) alongS¹(x⁰ x t).

The matrix

G²(x⁰ x t) := G¹(x⁰ x t) +f_x⁰(x⁰ x t)PQ¹(x⁰ x t) remains nonsingular on^GI.

Clearly, all matrices given here depend continuously on their arguments.

4

Now let tl ^{2 I} denote the current time and let x⁰_l xl be the given approximations of (Px)⁰(tl) x(tl), which cause the defect

l :=f(x⁰_l xl tl)⁶= 0: (2.3)

It should be emphasized once more that it does not matter at all howx⁰_l xlare computed.

Usually,x⁰_l represents a nite dierence, typicallyP_h¹l(xl^;xl^;1).

We aim at a new, improved approximation x^newl such that the critical part of the defect l vanishes. However, what is the critical part ofl and which component of the approxi- mations correspond to it?

For shortness, let us write : Q¹l :=Q¹(x⁰_l xl tl), P¹l :=I^;Q¹l, G²l :=G²(x⁰_l xl tl). We refer to Section 4 below for the numerical computation ofQ¹l.

By means of a careful error analysis of integration methods it is shown (e.g. Tischendorf (1996), Marz (1992)) that PQ¹lG^;12_l l represents the critical part of the defect. On the other hand, there is a close relationship between that part of the defect and the error of the PQ¹l-component of the solution.

Consider the ansatz

x^new_l =xl+PQ¹z (2.4)

PQ¹lG^;12_l f(x⁰_l x^new_l tl) = 0 (2.5)

z =PQ¹lz: (2.6)

In Section 3 it will be shown that (2.4) - (2.6) determinex^new_l z locally uniquely provided that the approximations we start with are suciently close. With the initial guessx^new(0)_l = xl z⁽⁰⁾ = 0 one Newton iteration for the system

x^new_l ^;xl^;PQ¹lz = 0

PQ¹lG^;12_l f(x⁰_l x^new_l tl) + (I^;PQ¹l)z = 0 yields

x_Nl :=x^new(1)_l =xl^;PQ¹lG^;12_l l (2.7) which is the projected defect correction formula we are looking for. In Section 4 it is described how to realize (2.7) numerically.

In the special case of linear DAEs

A(t)x⁰(t) +B(t)x(t)^;q(t) = 0 (2.8) we have, on the one hand,

A(tl)(Px)⁰(tl) +B(tl)x(tl)^;q(tl) = 0 but, on the other hand,

A(tl)x⁰_l+B(tl)xl^;q(tl) =l:

5

Taking into account the given relations A(tl) =G¹lP =G²lP¹lP

PQ¹lG^;12_l A(tl) =PQ¹lG^;12_l G²lP¹lP =PQ¹lP¹lP = 0 Q¹l =Q¹lG^;12_l B(tl) =Q¹lG^;12_l B(tl)P

we nd

PQ¹lxl =PQ¹lG^;12_l q(tl) +PQ¹lG^;12_l l

PQ¹lx(tl) =PQ¹lG^;12_l q(tl)

thusPQ¹lxl^;PQ¹lG^;12_l l=PQ¹lx(tl):Consequently, the defect correction formula (2.7) means for linear DAEs

x_Nl = (I^;PQ¹l)xl+PQ¹lxl^;PQ¹lG^;12_l l

= (I^;PQ¹l)xl+PQ¹lx(tl)

i.e., the PQ¹l-component of the approximation is replaced by the corresponding compo- nent of the true solution.

Note that in the linear case it holds that x^new_l =x_Nl:

Next we show the projected defect correction formula (2.7) to generalize (1.6). For that, we return to the Hessenberg DAE (1.1), (1.2) and formQ=diag(0 I),

PQ¹l =

Hl 0 0 0

!

PQ¹lG^;12_l =

0 Fl

0 0

!

:

Recall that Hl represents a projector. It projects onto im^@b_@x¹2(xl¹ xl² tl) along ker_@x^@b21(xl¹ tl).

In this case, formula (2.7) reads

x_Nl¹ =xl^1;Flb²(xl¹ tl) (2.9)

x_Nl² =xl² (2.10)

that is, again we have the projected defect correction formula (1.6) for Hessenberg form DAEs .

6

3 Error Analysis

Now we have to ask whether x^new_l is well-dened by the nonlinear equation system (2.4)- (2.6) and how it is related to the true solution valuex(tl). At the same time we try to give an estimate for the rst Newton-iterationx_Nl.

Let the DAE (1.2) be index-2 tractable in the neighbourhood^GI of the graph ; of the true solution as supposed in Section 2.

As above,tl^2I and (x⁰_l xl)^2G denote the current time and the given approximation to be improved, respectively. Q¹l Q²l etc. are the matrices used in Section 2.

Next introduce the auxiliary function

f~(x⁰ x) :=PQ¹lG^;12_l f(x⁰ x tl) (x⁰ x)^2G (3.1) that is continuously dierentiable due to our general assumptions agreed upon above. Of course, ~f depends on l, but we drop that index.

When investigating the solvability of DAEs with index 2, functions like ~f are supposed to be at least in C² (e.g. Tischendorf (1996)). From this point of view, it seems to be very natural that we now assume the Lipschitz-conditions

jf~_x⁰(x⁰ x)^;f~_x⁰(x⁰ x)^jL^1jx^;x^j+L^2jx^0;x^0j (3.2)

jf~_x⁰⁰(x⁰ x)^;f~_x⁰⁰(x⁰ x)^jL^3jx^;x^j+L^4jx^0;x^0j (3.3) (x⁰ x) (x⁰ x)^2G

to be satised.

It should be noted that our construction leads to the relation PQ¹l =PQ¹lG^;12_l f_x⁰(x⁰_l xl tl) = ~f_x⁰(x⁰_l xl)

which will be used later.

Furthermore, if we denote the critical part of the defect shortly by ~l, we have ~l :=PQ¹lG^;12_l l=PQ¹lG^;12_l f(x⁰_l xl tl) = ~f(x⁰_l xl):

In applications, the main interest is directed to special quasi-linear DAEs

a(x t)x⁰(t) +b(x(t) t) = 0: (3.4)

For those equations the auxiliary function ~f becomes invariant of its rst argument if the image space of the leading coecient matrixa(x t) does not rotate with changingx. This fact is gured out by the following lemma. Things become essentially simpler in that case.

7

Lemma 3.1

^{Let f(x0 x t)a(x t)x0+b(x t) and let im a(x t}) be independent of x. (i) Then it holds that

f~(x⁰ x)PQ¹lG^;12_l b(x tl):

(ii) The subspaces S(x⁰ x t) N¹(x⁰ x t) S¹(x⁰ x t) as well as the matrix G¹(x⁰ x t) do not depend on their rst argument x⁰ at all.

(iii) The projector Q¹(x⁰ x t) onto N¹(x⁰ x t) along S¹(x⁰ x t) is invariant of x⁰. (iv) The expressionQ¹G^;12 does not vary withx⁰, andQ¹G^;12 =Q¹(G¹+b⁰_xPQ¹)^;1 holds

true.

Proof:

Sinceim a(x t) does not vary withx, it follows that a⁰_x(x t)w²im a(x t) w²IR^m

holds true. Since Q is assumed to remain constant, the matrix function G¹ simplies to G¹ =a+b⁰xQ. Moreover, now we have

S¹(x⁰ x t) =^fz ²IR^m :b⁰_x(x t)Pz ²im(a(x t) +b⁰_x(x t)Q)^g:

By this, assertion (ii) becomes straight forward. Further, (ii) implies (iii) due to the uniqueness of that projection. Further, the matrix function ~G² :=G¹+b⁰_xPQ¹ is nonsin- gular and assertion (iv) follows immediately. To show (i), consider the relation

im a(xl tl)²kerPQ¹lG^;12_l

given by our construction. Namely, it holds thatQ¹lG^;12_l f_x⁰⁰(x⁰_l xl tl) =Q¹lG^;12_l G²lP¹lP = 0, thus Q¹lG^;12_l a(xl tl) = 0.

Because of im a(x tl) im a(xl tl) we nally obtain ~f(x⁰ x) PQ¹lG^;12_l (a(x tl)x⁰ + b(x tl))PQ¹lG^;12_l b(x tl): ²

Corollary 3.2

Lemma 3.1 impliesL² =L³ =L⁴ = 0 in the conditions (3.1), (3.2).

Now turn to the nonlinear system (2.4)-(2.6) that can equivalently be rewritten as

x^new_l =xl+z (3.5)

PQ¹lG^;12_l f(x⁰_l xl+z tl) + (I^;PQ¹l)z = 0: (3.6) Hence, we may solve equation (3.6) forz and then determine x^new_l by means of (3.5).

8

Theorem 3.3

Given suciently accurate approximationsx⁰_l xl of (Px)⁰(tl) x(tl): (i) Then there is a radius % > 0 such that equation (3.6) has a unique solution

z ²B(0 %).

(ii) For x^new_l :=xl+z the error estimation

jPQ¹l(x^new_l ^;x(tl))^jM^1jxl^;x(tl)^j2+M^2jx⁰_l^;(Px)⁰(tl)^j2 is valid with certain constants M¹ M².

Proof:

Introduce the function

F(z) := ~f(x⁰_l xl+z) + (I^;PQ¹l)z z ²B(0 %)

choosing% >0 small enough to realize L¹% <1 and ^f(x⁰_l xl+z) :^jz^j%^{g G}: Solving equation (3.6) means in fact to look for a zero ofF.

A priori, the functionF belongs toC¹. Compute F⁰(z) = ~f_x⁰(x⁰_l xl+z) + (I^;PQ¹l)

= PQ¹lG^;12_l (f_x⁰(x⁰_l xl+z tl)^;f_x⁰(x⁰_l xl tl)) +PQ¹l+ (I^;PQ¹l)

= I+r(z)

wherer(z) := ~f_x⁰(x⁰_l xl+z)^;f~_x⁰(x⁰_l xl tl) ^jr(z)^jL^1jz^j, z ²B(0 %): Further, the Lipschitz condition

jF⁰(z)^;F⁰(z)^jL^1jz^;z^j z z²B(0 %)

results immediately, and the Jacobian F⁰(z) remains nonsingular on B(0 %).

Next we show that Banach's Fixed Point Theorem applies to the map

Az :=z^;F(z) z²B(0 %)^\im PQ¹l =:^M(%): Forz z^2M(%) we derive

jAz^;Az^{j Z 1}

0

jI^;F⁰(sz + (1^;s)z)^jds^jz^;z^jL¹%^jz^;z^j and

Az = z^;f^~(x⁰_l xl+z)

= z^;( ~f(x⁰_l xl+z)^;f~(x⁰_l xl))^;~l

= z^;R01f^~_x⁰(x⁰_l xl+sz)ds z^;~l

= z^;R01ff^~_x⁰(x⁰_l xl+sz)^;f~_x⁰(x⁰_l xl)dsz^;PQ¹lz^;~l

thus

jAz^{j 1}₂L^1jz^j2+^j~l^j 1

2%+^j~l^j: 9

It becomes clear that, actually,^Amaps^M(%) into itself contractively if the given approx- imations are suciently accurate to satisfy the condition

j~l^j 1

2% (3.7)

i.e., assertion (i) is proved.

Next, consider the resulting error (cf. (3.5)) x^new_l ^;x(tl) =xl+z^;x(tl)

whereF(z) = 0 as well as xl and x(tl) are close to each other so that x(tl)²B(xl %).

Temporarily, we abbreviate x := x(tl) x⁰ := (Px)⁰(tl). Recall that x^new_l ^; xl = z PQ¹lz =z, F(x^new_l ^;xl) = 0 and derive

F(x^;xl) = F(x^;xl)^;F(x^new_l ^;xl)

= ^R01F⁰(sx + (1^;s)x^new_l ^;xl)ds(x ^;x^new_l ): Since the matrix

El :=^{Z 1}

0

F⁰(sx+ (1^;s)x^new_l ^;xl)ds (3.8)

satises the inequality ^jEl^;I^j L¹% < 1, it is invertible such that we obtain the error representation

x^;x^new_l =E_l^;1F(x^;xl): (3.9)

Due to its construction, El has the properties

PQ¹l(El^;I) =El^;I (I^;PQ¹l)El=I^;PQ¹l

which lead to

(I ^;PQ¹l)E_l^;1 = (I^;PQ¹l)

PQ¹lE_l^;1 = E_l^;1PQ¹l+E_l^;1(I^;El)(I^;PQ¹l): Together with (3.9) this implies

(I^;PQ¹l)(x^;x^new_l ) = (I^;PQ¹l)(x^;xl) (3.10) and

PQ¹l(x^;x^new_l ) =E_l^;1f^~(x⁰_l x) +E_l^;1(I^;El)(I^;PQ¹l)(x^;xl): (3.11) The rst relation (3.10) conrms that exactly thePQ¹l-component of the given approxi- mationxl is changed. Let us take a closer look at the two terms of the right-hand side in formula (3.11). We should consider in some detail thatEl depends on x^new_l .

10

We have

f~(x⁰l x) = ~f(x⁰l x)^;f~(x⁰ x) =

= ^R01f~_x⁰⁰(sx⁰_l+ (1^;s)x⁰ x)ds(x⁰_l^;x⁰)

= ^R01ff~_x⁰⁰(sx⁰_l+ (1^;s)x⁰ x)^;f~_x⁰⁰(x⁰_l xl)^gds(x⁰_l^;x⁰) hence

jf~(x⁰_l x)^j 1

2L^4jx⁰_l^;x^0j2+L^3jxl^;x^jjx⁰_l^;x^0j: (3.12) On the other hand, from

(El^;I)(I^;PQ¹l) =^{Z 1}

0

ff~_x⁰(x⁰_l sx+ (1^;s)x^new_l )^;f~_x⁰(x⁰ x)^g(I^;PQ¹l) we derive

j(El^;I)(I^;PQ¹l)(x^;xl)^jf12L^1jx^new_l ^;x^j+L^2jx⁰_l^;x^0jgj(I^;PQ¹l)(x^;xl)^j

f 1

2L^1j(I^;PQ¹l)(xl^;x)^j+L^2jx⁰_l^;x^0jgj(I^;PQ¹l)(x^;xl)^j+

+¹²L¹%^jPQ¹(x^new_l ^;x)^j (3.13)

where ^jI^;PQ¹l^j denotes a constant such that^j(I^;PQ¹l)(x^;xl)^j%. Using the bound^jE_l^;1j(1^;L¹%)^;1 and inserting (3.12), (3.13), into (3.11) we nd

jPQ¹l(x^new_l ^;x)^j(1^;L¹%)^;1f12L^4jx⁰_l^;x^0j2+L^3jxl^;x^jjx⁰_l^;x^0j+ +¹²L^1j(I^;PQ¹l)(xl^;x)^j2+L^2jx⁰_l^;x^0jj(I^;PQ¹l)(x^;xl)^jg+

+(1^;L¹%)^;1L¹%^12jPQ¹l(x^new_l ^;x)^j: ^(3.14) Finally, supposing further % to be small enough to satisfy %L¹(1 + ²)<1, we obtain the estimation oered in assertion (ii) with

M¹ := ¹² 1^;%L¹(1 + ²)^;1(L^1jI ^;PQ¹l^j2+L^2jI^;PQ¹l^j+L³) M² := ¹² 1^;%L¹(1 + ²)^;1(L²+L³+L⁴):

2

Corollary 3.4

In case of Lemma 3.1 assertion (ii) of Theorem 3.3 simplies to

jPQ¹l(x^new_l ^;x(tl))^j 1

2L¹1^;%L¹(1 + 2)

;1

j(I^;PQ¹l)(xl^;x)^j2: (3.15)

Proof:

The inequality (3.15) is a direct consequence of (3.14). ²

11

Remarks

1. Starting a Newton iteration for solvingF(z) = 0 with the initial guessz^{(0) 2}B(0 %), we nd

z^(1);z = z^(0);F⁰(z⁽⁰⁾)^;1F(z⁽⁰⁾)^;z

= F⁰(z⁽⁰⁾)^{;1Z 1}

0

fF⁰(z⁽⁰⁾)^;F⁰(sz⁽⁰⁾+ (1^;s)z)^gds(z^(0);z)

jz^(1);z^j(1^;L¹%)^{;1 12}L^1jz^(0);z^j2.

This shows B(0 %) to belong to the region where the Newton iterations converge quadratically.

For z⁽⁰⁾ = 0, F⁰(0) =I it results that z⁽¹⁾ =^;F(0) = ^;~l:

2. Summarizing the conditions for x⁰_l xl to be suciently accurate, we have standard restrictions L¹%(1 +²) <1 (which impliesL¹% < 1) and ^jx(tl)^;xl^j %, but also

j~l^{j 12}%.

In special cases the size of ~l is fully governed by xl (cf. Lemma 3.1) and the P-component of xl (Hessenberg form DAEs), respectively.

3. Of course, also the simple iterationsz^0]:= 0, z^j+1] :=z^j];F(z^j]) remain in^M(%) and converge to z.

Lemma 3.5

It holds that

~l =PQ¹l(xl^;x(tl)) +Rl (3.16)

with

jRl^{j 12}L¹+¹⁴(L²+L³)^jxl^;x(tl)^j2+ + ¹²L⁴ +¹⁴(L²+L³)^jx⁰_l^;(Px)⁰(tl)^j2

Proof:

Use once more the abbreviationsx :=x(tl),x⁰ := (Px)⁰(tl). Compute

~l = ~f(x⁰_l xl) = ~f(x⁰_l xl)^;f~(x⁰ x)

= ^R01nf~_x⁰⁰(sx⁰_l+ (1^;s)x⁰ sxl+ (1^;s)x)(x^0;x⁰_l) + ~f_x⁰(sx⁰_l+ (1^;s)x⁰ sxl+ (1^;s)x)(x ^;xl)^ods

= ^R01nf~_x⁰⁰(sx⁰_l+ (1^;s)x⁰ sxl+ (1^;s)x)^;f~_x⁰⁰(x⁰_l xl)^ods(x^0;x⁰_l) +^R01nf~_x⁰(sx⁰_l+ (1^;s)x⁰ sxl+ (1^;s)x)^;f~_x⁰(x⁰_l xl)^ods(x^;xl)

+PQ¹l(x^;xl)

= Rl+PQ¹l(x^;xl)

12

and estimate

jRl^{j 12}L^4jx^{0 ;}x⁰_l^j2 + ¹²(L³+L²)^jx^0;x⁰_l^jjx^;xl^j

+ ¹²L^1jx^;xl^j2:

As far as the errorx_Nl^;x(tl) is concerned, Lemma 3.5 indicates what it looks like. Hence,2

the next theorem is obvious.

Theorem 3.6

For the rst Newton iteration x_Nl given in (2.7) it holds that (I^;PQ¹l)(x_Nl ^;x(tl)) = (I^;PQ¹l)(xl^;x(tl))

and

jPQ¹l(x_Nl ^;x(tl))^jjRl^j:

Corollary 3.7

For quasi-linear DAEs (3.4) in case of Lemma 3.1 it holds that

jPQ¹l(x_Nl ^;x(tl))^j 1

2L^1jxl^;x(tl)^j2:

It should be mentioned that there are (Freude (1995)) certain similar convergence results as described in Theorem 3.3 concerning a more special class of index-2 DAEs, where formally the projector Q¹((Px)⁰(tl) x(tl) tl) and the corresponding matrices are used instead ofQ¹l =Q¹(x⁰_l xl tl) etc. in our context. Then, continuity arguments are applied to justify the practical use ofQ¹l.

Another straightforward generalization of (1.3)-(1.5) would use the new approximation x^new_l also for determining the back-projection. In this case, instead of (2.4)-(2.6), one has to solve the nonlinear system

x^new_l ^;xl^;PQ¹(x⁰_l x^new_l tl)z = 0 (3.17) (PQ¹G^;12 f)(x⁰_l x^new_l tl) + (I ^;PQ¹(x⁰_l x^new_l tl))z = 0 (3.18) which is even practically more expensive. One Newton step with the initial guess

x^new(0)_l :=xl z⁽⁰⁾ := 0 leads to

x^new(1)_l =xl+PQ¹lz⁽¹⁾

(PQ¹G^;12 f)⁰_x(x⁰_l xl tl)(x^new(1)_l ^;xl) + (I ^;PQ¹l)z⁽¹⁾ =^;~l: Considering

(PQ¹G^;12 f)⁰_x(x⁰_l xl tl) = PQ¹lG^;12_l f_x⁰(x⁰_l xl tl)+

+(PQ¹G^;12 )⁰_x(x⁰_l xl tl)f(x⁰_l xl tl) =

= PQ¹l+ (PQ¹G^;12 )⁰_x^jll

13

we know the Jacobian of (3.17), (3.18) to become nonsingular at the initial guess provided that the term (PQ¹G^;12 )⁰_x^jllis small enough, which can be ensured for suciently accurate approximations x⁰_l xl such thatl becomes small.

Hence, the rst Newton iteration yields

x^new(1)_l =xl+PQ¹lz⁽¹⁾ (3.19)

PQ¹l(x^new(1)_l ^;xl) + (I^;PQ¹l)z⁽¹⁾ =^;~l^;(PQ¹G^;12 )⁰_x^jllPQ¹lz⁽¹⁾ (3.20) that is,PQ¹lz⁽¹⁾ in (3.19) is now the solution of the linear equation

PQ¹lz⁽¹⁾ =^;~l^;(PQ¹G^;12 )⁰_x^jllPQ¹lz⁽¹⁾: (3.21) From this point of view, the projected defect correction formula (2.7) represents an in- complete realization of (3.19), (3.20) that corresponds to a single simple iteration step in the linear equation (3.21) to be solved forz⁽¹⁾.

Finally, emphasize once more that the matrix (PQ¹G^;12 )⁰_xl seems not to be available in practice in general.

4 Computing the projector

^Q1l

numerically

To realize the projected defect correction (2.4)-(2.6) or (2.7), the projector valueQ¹l :=

Q¹(x⁰_l xl tl), that is, the projector onto the nullspace N¹(x⁰_l xl tl) along the associated subspace S¹(x⁰_l xl tl), is needed. Both subspaces are determined by the matrix pair

fCl Dl^g via

Cl := f_x⁰⁰(x⁰_l xl tl) +f_x⁰(x⁰_l xl tl)Q (4.1) Dl := f_x⁰(x⁰_l xl tl)P

N¹(x⁰_l xl tl) = kerCl

S¹(x⁰_l xl tl) = ^fz ²IR^m :Dlz ²im Cl^g: (4.2) The projectorQ used to formCl Dl is often known a priori by the DAE structure. Typi- cally, for semi-explicit DAEs, we simply haveQ=diag(0 I). IfQis not given in advance, it can be provided easily, say by a Householder decomposition. Note that any projector ontoN can be used for Q, i.e., there is no need to have the orthoprojector.

Looking at (4.2), (4.2) we know the projectorQ¹l to be nothing else but the canonical pro- jector of the matrix pair ^fCl Dl^g (cf. Griepentrog and Marz (1986), page 201). Namely, because of the index-2 condition (2.2) for the DAE (2.1), the matrix pair ^fCl Dl^g has index 1. Hence, the canonical projector of^fCl Dl^gis represented by

Vl(Cl+DlVl)^;1Dl (4.3)

whereVl denotes any projector onto the nullspace ofCl. Recall that the matrixCl+DlVl

remains nonsingular due to the index-1 property of the pair^fCl Dl^g. 14