Asymptotic properties of regularized di erential-algebraic equations
Michael Hanke Antonio R. Rodrguez S.
1 Introduction
In the present paper we consider asymptotic properties of a regularization method for linear fully implicit di erential-algebraic equations (dae's)
A(t)x0+B(t)x(t) = 0 t 2t0 1): (1.1) We will assume that A(t) is singular for all t2t0 1). Such problems arise naturally in a number of applications, e.g. electrical networks, constraint mechanical systems of rigid bodies, chemical reaction kinetics at least as linearizations of nonlinear problems. This is why great interest has been devoted to the analysis, geometry and numerical treatment of dae's in recent years.
Nowadays, it is well-known that dae's (1.1) incorporate a great deal of new features compared with explicit ordinary di erential equations (i.e. equations solved for x0). A rough criterion for distinguishing between di erent classes of dae's is given by the notion of an index of a dae. Although there are a number of di erent notions (with varying aims and applicability) in common use, their common aim is to classify how far a given dae di ers from an explicit ordinary di erential equation. Ordinary dae's (i.e. A(t) is nonsingular for all t 2 0 1)) are characterized by the index 0. The index 1 describes the "most simple" class of dae's. Higher index equations are those with an index greater than 1. Analytically, the algebraic relations contained in (1.1) cause the higher index equations to include di erentiation problems such that nite di erence methods become unstable. This instability may become very dangerous if the nullspace of A(t) varies with
t (9]). On the other hand, a number of numerical methods is available for the solution of index 1 dae's. So a possible way to solve higher index problems is the index reduction.
In the present paper we consider a regularization method. Eq. (1.1) is perturbed by a small parameter" such that the resulting system has index 1 and, for "!0, the solution of the regularized systems tends to that of (1.1). We consider the approach of 7] whose convergence properties are analyzed in 4, 5]. Note that this approach is closely related to other parametrizations, among them Baumgarte stabilization (cf. 1]).
We are interested in stability properties of the regularization approach. Namely, we will show that if (1.1) is asymptotically stable (in a sense given below), then the regularized systems are so, too. Similar results for autonomous quasilinear systems are given in 10]
relating on stability criteria proved in 2] and 8]. We emphasize that our results are valid if the nullspace of A(t) depends on t.
The paper is organized as follows. In Section 2 we introduce the solution representation of linear index 1 and index 2 dae's. Motivated by the notion of exponential asymptotic
1
stability 3] we generalize this notion to the dae case. Section 3 deals with the regulariza- tion of linear index 2 dae's and some of its properties. Finally, in Section 4, we present our main result concernig the exponential asymptotic stability of the regularization.
2 Exponential asymptotic stability for linear dae's
In a rst step consider the explicit ordinary di erential equation
x
0=B(t)x t2t0 1) (2.1)
with continuous coecients B 2C(t0 1) L(IRm)),x(t)2IRm.
De nition 1
(cf. 3, p. 84])The trivial solution of (2.1) is called exponentially asymptotically stable (eas) if there are constants K > 0 such that, for all t t0, x0 2 IRm, the solution of the initial value problem
x
0+B(t)x= 0 t2t 1)
x(t) =x0 fulls the estimate
jx(t)jKjx0je; (t;t) t0 t t<1:
Remark:
(i) If the trivial solution of (2.1) is eas, then it is asymptotically stable (in Lyapunov's sense).
(ii) If B(t)B is a constant matrix function, then it is eas if and only if the real parts of the eigenvalues of B are strictly negative.
It is convenient to generalize Denition 1 slightly.
De nition 2
For every t 2 t0 1), let V(t) denote a subspace of IRm. The trivial solution of (2.1) is called exponentially asymptotically stable with respect to V if there are constants K >0 such that, for all t t0,x0 2V(t), the solution of the initial value problemx
0+B(t)x= 0 t2t 1)
x(t) =x0 fulls the estimate
jx(t)jKjx0je; (t;t) t0 t t<1:
In contrast to Denition 1 the set of admissible initial values is restricted. Denition 2 seems to be formaly. However it is useful if the subspaces V(t) are invariant solution spaces, i.e., if x(t) is a solution of (2.1) on t 1) with x(t) 2 V(t), then x(t) 2 V(t) for allt 2t 1).
2
Consider now the linear equation
A(t)x0+B(t)x= 0 t 2t0 1) (2.2) with continuous coecients. Assume that the nullspaceN(t) ofA(t) is smooth, i.e. there exists a continuously di erentiable matrix function Q 2 C1(t0 1) L(IRm)) such that
Q(t) is a projection onto N(t). Note that the rank of A(t) is constant then. The trivial case N(t)f0g is equivalent to (2.1) and shall be excluded. Furtheron, let P =I;Q. For simplicity we will drop the argumentt if no confusion can arise.
The nullspaceN of the leading coecient matrixA(t) determines which kind of functions we should accept for solutions of (1.1). Namely,AQ0 impliesAx0=APx0=A((Px)0;
P 0
x). Thus, (2.2) may be rewritten
A(Px)0+ (B;AP0)x= 0 t 2t0 1): (2.3) Therefore, we are looking for solutions of (2.2) in the function space
C
N1t0 1) :=fy2Ct0 1)jPy 2C1t0 1)g:
More precisely, a function x : t0 1) ! IRm is called a solution of (2.2) if it belongs to
C
N1t0 1) and fulls (2.1). Let
B
0 :=B;AP0 A1 :=A+B0Q:
Ifx is a solution of (2.2), then, for all t2t0 1), x(t) belongs to
S(t) :=fz 2IRm jB0(t)z2R(A(t))g
where R(A(t)) denotes the range of A(t). (2.2) is called transferable (or index 1 dae) if, for all t2t0 1),
N(t)S(t) =IRm: (2.4)
(2.4) is equivalent to the condition that A1(t) is nonsingular (cf. 2, Theorem A.13]).
Obviously, (2.3) is equivalent to
A
1
fP(Px)0+Qxg+B0Px= 0: (2.5) If (2.2) has index 1, we obtain by multiplication of (2.5) byQA;11 andPA;11 , respectively, the equivalent system
(Px)0;P0Px+PA;11 B0Px = 0
Qx+QA;11 B0Px = 0:
Lemma 1
If (2.2) has index 1, then it is equivalent to the systemu
0+ (PA;11 B0;P0)u = 0
v+QA;11 B0u = 0 t2t0 1) (2.6) where u = Px and v = Qx. Moreover, if u0 2 R(P(t)) for some t 2 t0 1), then the solution u of the initial value problem u0+ (PA;11 B0 ;P0)u = 0, u(t) = u0 fulls
u(t)2R(P(t)), t 2t0 1).
3
The proof of Lemma 1 is obvious. In order to show the second assertion, simply multiply the di erential equation by Q leading to (Qu)0;Q0Qu= 0.
Especially, (2.6) implies x = Px +Qx = u ; QA;11 B0u = (I ; QA;11 B0)Pu. Note that Qs(t) := QA;11 B0(t) is the projection of IRm onto N(t) along S(t). Once the Px- component of a solution is found, the nullspace components are given by a simple as- signment. Therefore, initial conditions can only be given for Px(t). This suggests the following denition.
De nition 3
Let (2.2) be an index 1 dae. The trivial solution of (2.2) is called expo- nentially asymptotically stable if there are constants K > 0 such that, for all t t0,x 0
2IRm, the solution of the initial value problem
A(t)x0+B(t)x= 0 t2t 1)
P(t)(x(t);x0) = 0 fulls the estimate
jx(t)jKjP(t)x0je; (t;t) t0 t t<1:
Remark:
In 2, p. 74] asymptotic stability in the sense of Lyapunov is dened for general nonlinear transferable dae's. If the trivial solution of (2.2) is eas, then it is asymptotically stable in the sense of the latter denition.Theorem 1
Let Qs =QA;11 B0 be bounded on t0 1). Then the trivial solution of (2.2) is eas if and only if the trivial solution of u0+ (PA;11 B0;P0)u = 0, t 2 t0 1), is eas with respect to R(P(t)).The proof follows simply from (2.6) and the representationx= (I;QA;11 B0)Pu.
Remark:
In 2, p. 78] a notion of contractivity for nonlinear transferable dae's is dened.Lemma 1.2.44 of that monograph shows that contractivity implies exponential asymptotic stability.
In contrast to (2.4), higher index dae's are characterized by nontrivial intersectionsN(t)\
S(t). Equivalently,A1 is singular. LetQ1(t) denote a projection onto the kernel ofA1(t),
B
1 =B0P A2 =A1+B1Q1:
The dae (2.2) is called tractable with index 2 if, for all t 2t0 1), A1(t) is singular but
A
2(t) is nonsingular. In this case (2.5) is equivalent to
A
2P1fP(Px)0+Qxg+Q1x] +B1P1x = 0: (2.7) Due to Lemma A.13 of 2]Q1A;12 B1(t) is a projection onto N(A1(t)) again. Assume that we have chosen this special projection from the beginning, i.e. Q1 Q1A;12 B1. Especially,
Q
1
Q = 0 holds such that Q, PQ1, PP1 are projections again and the pairwise product of two of them vanishes. Then, x can be decomposed like x = Qx +PQ1x+PP1x. Multiplying now (2.7) byPP1A;12 ,PQ1A;12 and QP1A;12 , respectively, we obtain
(PP1x)0+ (PP1A;12 B1;(PP1)0)PP1x = 0
PQ
1
x = 0
Qx+ ((QQ1)0+QP1A;12 B1)PP1x = 0
provided that Q1 is continuously di erentiable. The counterpart of Lemma 1 is now 4
Lemma 2
Let (2.2) be tractable with index 2, where Q1 Q1A;12 B1 is continuously dierentiable. Then (2.2) is equivalent to the systemz
0+ (PP1A;12 B1;(PP1)0)z = 0
y = 0 (2.8)
v = ;((QQ1)0+QP1A;12 B1)z
where z = PP1x, y = PQ1x and v = Qx. Moreover, if z0 2 R(PP1(t)) for some
t2t0 1), then the solution z of the initial value problemz0+(PP1A;12 B1;(PP1)0)z= 0,
z(t) =z0 fulls z(t)2R(PP1(t)), t 2t0 1).
Especially we have x = z +y +v = (I ;(QQ1)0 ;QP1A;12 B1)PP1z, where z solves an explicit ode with initial conditions z(t) 2 R(PP1(t)). Note again that (t) = (I ; (QQ1)0(t);QP1A;12 B1(t))PP1(t) is also a projector withN((t)) =N(PP1(t)). Similarly to the index 1 case, initial conditions can only be given for PP1(t). We are led to the following denition.
De nition 4
Let (2.2) be tractable with index 2, where Q1 =Q1A;12 B1 is continuously di erentiable. The trivial solution of (2.2) is called exponentially asymptotically stable if there are constants K >0 such that, for all tt0,x0 2IRm, the solution of the initial value problemA(t)x0+B(t)x= 0 t2t 1)
PP
1(t)(x(t);x0) = 0 fulls the estimate
jx(t)j KjPP1(t)x0je; (t;t) t0 tt<1:
Lemma 2 and Denition 4 yield immediately the following theorem.
Theorem 2
Let the assumptions of Lemma 2 be fullled. Moreover, let (QQ1)0 andQP
1 A
;1
2 B
1 be bounded. Then the trivial solution of (2.2) is eas if and only if the trivial solution of
z
0+ (PP1A;12 B1;(PP1)0)z = 0 t 2t0 1) is eas with respect to R(PP1(t)).
Remark:
(i) Again, if the trivial solution of (2.2) is eas, it is stable in the sense of Lyapunov.
(ii) Statements concerning Lyapunov stability of autonomous quasilinear index 2 and 3 systems can be found in 8].
(iii) Under the assumptions of Theorem 2, contractivity as dened in 6] is a sucient condition for eas of the trivial solution.
5
3 Regularization of index 2 dae's
A well-known method for approximating higher index dae's is the index reduction method.
We will consider a method based on regularization. More precisely, the index 2 dae (2.3) is replaced by the perturbed dae 4]
(A+"B0P)(Px)0+B0x= 0: (3.1) It is easy to see that for suciently small " 6= 0, (3.1) is an index 1 dae if (2.2) is tractable with index 2. For that denote A" = A +"B0P. For the nullspace we have
N(A(t)) = N(A"(t)) for all t 2 t0 1). Obviously, N(A(t)) N(A"(t)). If, for some
t 2 t0 1), z 2 N(A"(t)), then 0 = A"z = (A2+"B0P ;B0Q;B0PQ1)z. Multiplying
this equation byQ1A;12 yields"Q1z = 0. (Note thatQ1 =Q1A;12 B0P,A;12 B0Q=Q.) On the other hand, by multiplication by PP1A;12 one obtains PP1z+"PP1A;12 B0PP1z = 0.
If PP1A;12 B0 is uniformly bounded on t0 1), PP1z = 0 for all suciently small ". But this is equivalent to z 2 N(A(t)) because of PQ1z = 0. The relevant matrix for proving the transferability of (3.1) is now
A
1"=A"+B0Q=A+"B0P +B0Q: (3.2)
The nonsingularity of A1" can be shown similarly. Assume A1"z = y to hold, for xed
t 2 t0 1). Since A1" = A2;B0PQ1 +"B0P, we obtain by multiplication by Q1A;12 ,
PP
1 A
;1
2 and QA;12 , respectively,
"Q
1
z =Q1A;12 y
(I+"PP1A;12 B0)PP1z+PP1B0PQ1z =PP1A;12 y
Qz;QQ
1
z+"QA;12 B0P(P1+Q1)z=QA;12 y:
For given y 2 IRm, the rst equation gives Q1z, the second can be solved for PP1z provided that " is suciently small and PP1A;12 B0 is uniformly bounded. The nullspace component Qz is given by the third equation. This proves the desired nonsingularity.
Note that kA1"(t);1k=O(";1) since Q1z = 1"Q1A;12 y.
The convergence behaviour of solutions of initial value problems for (3.1) towards that of (2.2) can be characterized by asymptotic expansions on compact intervals. To be more precise, consider (2.2) on a compact interval t0 T] together with the initial conditions
PP
1(t0)(x(t0);x0) = 0: (3.3) Since (3.1) has index 1, we need additional initial conditions for PQ1x(t0) in order to specify a unique solution. In view of (2.8) it is reasonable to choose PQ1x(t0) = 0. For suciently small">0, letx" denote the solution of (3.1) subject to the initial condition (3.3) and
PQ
1
x(t0) = 0 (3.4)
on t0 T]. In 4] it is shown that, under suitable smoothness assumptions, an asymptotic expansion of the form
x"(t) =XN
j=0(xj(t) + xj())"j+O("N+1) (3.5) 6
with = t=", jxj()j Ce; holds true. x0 is the solution of (2.2) subject to (3.3).
If PQ1x(t0) is not chosen consistently (i.e. (3.4)), an additional term x;1()";1 with
Px;1 0 appears. The form of (3.5) leads to the conjecture that, if the index 2 problem (2.2) is asymptotically stable in some sense, then this should be true for the regularized system, too. Unfortunately in general such a statement does not hold. In the next section we will give sucient conditions for preserving exponential asymptotic stability. Since a number of regularizations introduced by other authors are closely related to (3.1), (1]), for them similar results should be expected.
It is convenient to decompose (3.1) as (2.8) (cf. 4]).
Lemma 3
Let the assumptions of Lemma 2 be fullled. Moreover, let C1 :=PP1A;12 B1 be continuously dierentiable andC2 :=QP1A;12 B1. Then(3.1) is equivalent to the system(I+"C1)z0+ (C1;(PP1)0)z;((PP1)0+"C10)PQ1y = 0
"y 0
;"(PQ1)0z+ (I;"(PQ1)0)y= 0 (3.6)
v=;C2z;"C2(z0+y0) +QQ1(z0+y0)
where z =PP1x, y=PQ1x, v =Qx. Moreover, if z0 2R(PP1(t)) for some t 2t0 1), then the solution of the initial value problem (I +"C1)z0+ (C1 ;(PP1)0)z ;((PP1)0+
"C 0
1)PQ1y = 0, z(t) = z0 fulls z(t) 2 R(PP1(t)) for t 2 t0 1) and any continuous function y. A similar assertion holds for y as a solution of the second equation.
4 Exponential asymptotic stability of regularized lin- ear dae's
The aim of this section is to show that exponential asymptotic stability of (2.2) carries over to (3.1) under some smoothness and boundedness assumptions.
Lemma 4
LetM(t)2C1t0 1) such that M andM0 are bounded. Then, for suciently small ", it holds:(i) X"(t) :=I +"M(t) is nonsingular.
(ii) X";1(t) is bounded on t0 1) and the bound does not depend on ".
(iii) kdtdX";1(t)k"K uniformly on t0 1).
Proof:
(i) and (ii) are simple consequences of Banach's theorem. (iii) follows from ddtX";1(t) =;X";1(t) ddtX"(t)X";1(t) =;"X"(t);1M0(t)X";1(t):Lemma 5
Consider the systemx
0+A(t ")x= 0 t 2t0 1) (4.1) for 0<""0 and the perturbed system
x
0+A(t ")x=F(t ")x t2t0 1): (4.2) Let A F : t0 1)0 "0]!L(IRm) be continuous matrix functions. Suppose that
7
(i) (4.1) is eas with respect to the set of subspaces V with constants ("), K("), (ii) R(F(t "))V(t) for every t2t0 1) and every " 2(0 "0],
(iii) kF(t ")k"K^ uniformly in t0 1) for all "2(0 "0], and (iv) (")>"K(") ^K, "2(0 "0].
Then (4.2) is eas with respect to V, and it holds that
jx(t)jK(")e;^ (")(t;t)jx(t)j t0 tt<1 "2(0 "0]
for the solution of (4.2) subject to the initial condition x(t) = x0 2 V(t). Here, ^(") =
(");"K(") ^K.
Proof:
Let (t s ") denote the fundamental solution of (4.1) subject to the initial con- dition (s s ") =I. Due to assumption (i) we havej (t s ")cjK(")jcje; (")(t;s) t0 st<1 (4.3) for c 2V(s). The solution of (4.2) subject to the initial condition x(t) = x0 is uniquely determined as a solution of the integral equation
x(t) = (t t ")x0+Zt
t (t s ")F(s ")x(s)ds tt:
Since F(s ")x(s)2V(s) by assumption (ii), (4.3) yields, for x0 2V(t),
jx(t)j K(")jx0je; (")(t;t)+Zt
t
K(")kF(s ")kjx(s)je; (")(t;s)ds
K(")jx0je; (")(t;t)+"K(") ^KZt
t
e
; (")(t;s)
jx(s)jds: (4.4) The latter inequality will be multiplied by e (")t. With the notation,
y(t) =Zt
t e
(")s
jx(s)jds we obtain
y 0
K(")jx0je (")t+"K(") ^Ky or, equivalently,
y 0
;"K(") ^KyK(")jx0je (")t: Multiplication by e;"K(")Kt^ and integration yields
y(t) jx0j
"K^e (")t(e;"K(")K^(t;t);1): Introduce this estimate in (4.4) to obtain
jx(t)jK(")jx0je;( (");"K(")K^)(t;t): This gives the assertion.
We are now ready to prove our main result.
8
