Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Numerical Methods for Stiff Index-3 DAEs

Numerical Methods for Stiff Index-3 DAEs Mathematical and Computer Modelling of Dynamical Systems 1387-3954/01/0702-239$16.00 2001, Vol. 7, No. 2, pp. 239±262 Swets & Zeitlinger I. HIGUERAS ABSTRACT Space semidiscretization of PDAEs, i.e. coupled systems of PDEs and algebraic equations, give raise to stiff DAEs and thus the standard theory of numerical methods for DAEs is not valid. As the study of numerical methods for stiff ODEs is done in terms of logarithmic norms, it seems natural to use also logarithmic norms for stiff DAEs. In this paper we show how the standard conditions imposed on the PDAE and the semidiscretized problem are formally the same if they are expressed in terms of logarithmic norms. To study the mathematical problem and their num- erical approximations, this link between the standard conditions and logarithmic norms allow us to use for stiff DAEs techniques similar to the ones used for stiff ODEs. The analysis is done for problems which appear in the context of elastic multibody systems, but once the tools, i.e., logarithmic norms, are developed, they can also be used for the analysis of other PDAEs / DAEs. Keywords: DAEs, elastic multibody systems, logarithmic norms, Runge±Kutta method, stiff. 1 INTRODUCTION Nowadays there is an increasing interest in the analysis and ef®cient numerical solution of coupled systems of partial differential equations (PDEs) and algebraic equations. This kind of problems appear for example in models of constrained mechanical systems including rigid and elastic bodies [11]. In some ®elds, these coupled systems are known as partial differential-algebraic systems (PDAEs). The aim of this paper is to analyze and obtain valid bounds for the numer- ical solution of PDAEs / DAEs which appear in the context of mechanical systems. We will deal with problems of the form 00 0 u ‡ Au‡ B  ˆ l …1† Bu ˆ m …2† Departamento de Matema Âtica e Informa Âtica, Universidad Pu  blica de Navarra, 31006 Pamplona, Spain. E-mail: higueras@unavarra.es 240 I. HIGUERAS 0 0 with A and B time independent linear maps, A : V ! V , B : V ! Q , V and Q 0 0 0 are Hilbert spaces, V and Q are their dual spaces respectively, and B denotes the dual map of B. The map A is such that there exists a constant > 0 with hAu; ui ;kuk for all u 2 Ker B ; …3† and B such that there exists a constant > 0 with kB k kk for all  2 Q : …4† n m 0 t For the ®nite dimensional case, we will take V ˆ R , Q ˆ R , and B ˆ B . Recall that we have in mind not only DAEs in the standard case but also DAEs in the general operator sense. The reason to consider only these simple systems is that they cover both PDAEs formulated as operator DAEs and DAEs obtained from a PDAE by a semidiscretization process in space. As these problems are stiff or may contain unbounded operators, the standard theory of numerical methods for DAEs [6, 9] is not valid and a new analysis must be done. Taking into account that in the theory of numerical methods for stiff ODEs the tool used is the logarithmic norm of a map (more precisely, the l.u.b. logarithmic norm), it seems natural to use logarithmic norms also for stiff DAEs. This paper shows how the standard conditions imposed on the operator formulation of the PDAE and the semidiscretized problem are formally the same if they are expressed in terms of logarithmic norms. To study the mathematical problem and the numerical approximations, this link between the standard conditions and logarithmic norms allows us to use for stiff DAEs similar techniques to the ones used for stiff ODEs. In Section 1, we give the PDAE formulation as operational DAEs and the semidiscretization process, and show that (1)±(2) together with the conditions (3)±(4) cover both problems. The section ends with a brief outline of numerical methods for stiff ODEs. In Section 2, devoted to logarithmic norms, we summarize the main de®nitions, extensions and results on this topic. The origin for this paper are problems which appear in the context of elastic multibody systems (1)±(2), but we also have in mind more general DAEs. This is the reason to extend the concept of g.l.b. logarithmic norm of a matrix to a matrix pencil following the lines of the one done in [7] for the l.u.b. logarithmic norm of a matrix. Once we have constructed the tools, we can use them to analyze the system (1)±(2) and its numerical solution with stif¯y accurate Runge±Kutta methods. This is done in Section 4. Finally, in Section 5 some conclusions and forthcoming work is pointed out. STIFF INDEX-3 DAEs 241 2 MATHEMATICAL MODELS AND NUMERICAL METHODS BACKGROUND Mathematical modeling may lead to coupled systems of partial differential equations and algebraic equations. In [11], in the context of elastic multibody systems, the following PDAE in weak formulation is studied, ``Find …u…; t†;…; t†† 2 V  Q such that hu ; vi‡ a…u; v†‡ b…v;†ˆhl; vi for all v 2 V …5† b…u;†ˆhm;i for all  2 Q '' ; …6† where V and Q are two Hilbert spaces, and a : V  V ! R, and b : V  Q ! R are two continuous bilinear forms. If we denote V ˆfv 2 V : b…v;†ˆ 0 for all  2 Qg ; it is assumed that the bilinear form a…;† is V -elliptic, i.e., there exists a constant > 0 such that a…u; u† kuk for all u 2 V : Moreover, the bilinear form b…;† is assumed continuous, i.e., there exists a constant  such that ka…u; v†k  kukkvk for all u 2 V : A standard condition for the bilinear form b…;† is the LBB-condition, i.e., there exists a constant > 0 such that b…v;† inf sup  : 2Q kvk kk v2V V Q 0 0 As it is pointed out in [11], we can consider the operators A : V ! V , with V the dual space, de®ned as hAu; viˆ a…u; v†, B : V ! Q , de®ned as 0 0 hBv;iˆ b…v;† and its dual operator B : Q ! V , de®ned as hB ; viˆ b…v;†. In this way, (5)±(6) is 00 0 u ‡ Au‡ B  ˆ l …7† Bu ˆ m …8† In terms of the map B, we have that V ˆ Ker B. The V -ellipticity of the 0 0 bilinear form a…;† can be written as hAu; ui kuk for all u 2 Ker B 242 I. HIGUERAS with > 0, the continuity of the bilinear form a…;† as kAuk kuk for all u 2 V ; and the LBB-condition as kB k kk for all  2 Q with > 0. Following the PDE approach, a method to solve numerically PDAEs consists in a semidiscretization in space followed by a time integration. In order to semidiscretize in space the equations (5)±(6), ®nite-dimensional subspaces V  V and Q  Q are considered, and solutions u …; t†2 V and  …; t†2 Q of the form q n X X u …x; t†ˆ v …x†q…t†; …x; t†ˆ …x†…t† i i  i i iˆ1 iˆ1 are constructed. If given a PDE the semidiscretization process gives an ODE system, for PDAEs, the problem obtained is a DAE system. The time dependent functions q…t†, …t† must solve the DAE i j 00 t M q …t†‡ A q…t†‡ B …t†ˆ l …t†…9† B q…t†ˆ m …t†…10† The mass matrix M and the matrix A are symmetric positive and semide®nite positive respectively. In (9)±(10) the matrices are de®ned in terms of the inner products in the Hilbert space and the bilinear forms, …M † ˆhv ; vi…A † ˆ a…v ; v†…B † ˆ b…v ;†: i j  i j  j i ij ij ij A compatibility condition is imposed on the discretization, 2 Q with b…v;†ˆ 0 for all v 2 V )  ˆ 0 which in terms of B is kB k kk with > 0 : With this condition the matrix B has full row rank and therefore the linear DAE (9)±(10) has index 3. To analyze convergence, the following subspace is considered W ˆf v 2 V : b…v;†ˆ 0 for all  2 Q g ; STIFF INDEX-3 DAEs 243 and it is imposed that the bilinear form a…;† is W -elliptic. Under these conditions, the matrix A satis®es hA q; qi kqk for all q 2 Ker B with > 0. The next step is to integrate (9)±(10) with some DAE method. We can expect that the nature of these DAEs be similar to the ODEs which result form the PDEs case: stiff problems. Therefore, the standard bounds which ensure certain order of convergence for index-3 DAEs [6, 9] are not valid anymore. Some authors [10] have studied the numerical solution of stiff mechanical systems using some model problems in which the stiffness parameter appears explicitly. As it was expected from the knowledge of stiff ODEs, depending on the stiffness parameter, an order reduction may appear. In the theory of numerical methods for stiff ODE problems, in order to analyze the effective order, the concept of B-convergence is introduced [2]. Basically the order of B-convergence is the order of the method independently of the stiffness of the problem. Other concepts used in this framework are the B-consistency and C-stability, which again are consistency and stability independently of the stiffness. Now the result is that B-consistency and C- stability imply B-convergence [2, Th. 7.2.5], and conditions to obtain C- stability and B- consistency must be studied. The class of ODEs for which this theory is developed is y ˆ f…t; y† such that h f…t; y †ÿ f…t; y †; y ÿ y i ky ÿ y k ; …11† 1 2 1 2 1 2 where h;i is an inner product norm in R . The class F is the class of problems satisfying (11) and such that their solutions are smooth, in the sense that the derivatives are bounded by constants independent of any quantity which may be in¯uenced by the stiffness [2, (7.2.2)]. The constant  is also known in the literature as the logarithmic norm of f , or the l.u.b. Dalhquist constant [2, 13]. When implicit Runge±Kutta methods are used, some nonlinear systems have to be used. In order to obtain existence and uniqueness of solution for these systems, the case < 0, is extremely favorable. Theorem 2.1 [2, Cor. 5.3.13] If a Runge±Kutta method is algebraically stable and irreducible, then the system of algebraic equations has unique solution for any problem from F with < 0. The important point in this result is that there is not any stepsize restriction. If > 0, a similar result is obtained with a stepsize restriction. 244 I. HIGUERAS Given a Runge±Kutta method with matrix coef®cients A, and in [2, (5.1.21)] the function hA ; i …A† ˆ min 6ˆ0 h;i is de®ned, where h;i denotes the inner product ~ ~ ~ h; i :ˆ d   ;  2 R ; d > 0 ; j ˆ 1; .. . ; s : j j j j jˆ1 In terms of this function, the following result is proved. The order of convergence is obtained from Theorem 7.3.3 in the same reference. Theorem 2.1 [2, Th. 7.4.4] Let a Runge±Kutta method…A; b † with positive weights be algebraically stable. Suppose that a positive diagonal matrix D exists such that …A† > 0. If q is the order stage of the method, then the method is B-convergent with order at least q on F . Again, the case < 0 is favorable in the sense that all the stability bounds involved to prove this result hold without any stepsize restriction. As a corollary, we obtain that the Gauss (order stage s), RadauIA (order stage sÿ 1), RadauIIA (order stage s), and the two stage Lobatto IIIC (order stage sÿ 1) methods are B-convergent. B-convergence of some other diagonally implicit and singly implicit methods are also proved in [2, Cor. 7.4.6, 7.4.7]. The numerical experiments shows that the lower bound for the order of convergence, i.e., the order stage of the method, can be reached. Recall that these proofs are done in terms of logarithmic norms and inner product norms, and they only use the standard properties of these concepts. 3 LOGARITHMIC NORMS The concept of logarithmic norm for a matrix was introduced in 1958 by Dahlquist and Lozinskij as a tool to study the growth of solutions to ODEs and the error growth in discretization methods for their approximate solution. For details see [2]. For a matrix A, the logarithmic norm (more precisely l.u.b. logarithmic norm) is de®ned by L‰I ‡  AŠÿ 1 M‰AŠˆ lim : …12† !0 STIFF INDEX-3 DAEs 245 The norm here is generally assumed to be an l.u.b. norm (or an l.u.b. Lipschitz constant), kAxk L‰AŠˆ max : …13† x6ˆ0 kxk Properties of norms and logarithmic norms can be found in [2, 1]. In particular, if the used norm is an inner product norm, we have the equivalent de®nition hAu; ui M‰AŠˆ max : …14† u6ˆ0 hu; ui In a similar way we can de®ne the g.l.b. logarithmic norms as l‰I ‡  AŠÿ 1 m‰AŠˆ lim ; …15† !0 where l‰Š is the g.l.b. Lipschitz constant kAxk l‰AŠˆ min : …16† x6ˆ0 kxk It can be proved that for the euclidean norm, p l ‰AŠˆ minf  :  eigenvalue of A Ag 2 i i ˆ minf :  singular value of Ag : i i Observe that from the de®nition l‰AŠ minfjj : det…I ÿ A†ˆ 0g : Moreover it can be proved that 0 if A singular l‰AŠˆ …17† ÿ1 ÿ1 L‰A Š if A regular and that m‰AŠˆÿM‰ÿAŠ. Thus if the norm is an inner product norm, hAu; ui m‰AŠˆ min : …18† u6ˆ0 hu; ui With the l.u.b. and g.l.b. logarithmic norms of a matrix, given the ODE y ˆ A…t†y, one can derive bounds for the solution and therefore give suf®cient 246 I. HIGUERAS conditions for stability. Some well-known theorems in stability theory of ODEs are easily proven with the help of logarithmic norms [15]. Remark 3.1 In the literature, usually these concepts are de®ned for square matrices but they can also be de®ned for rectangular matrices. & Remark 3.2 Observe that condition (4) can be expressed in terms of the g.l.b. Lipschitz constant l‰BŠˆ > 0. Hence we have t 2 l‰BBŠˆ > 0 ; and thus BB is a regular matrix. & Logarithmic norms for nonlinear maps are considered in [2, 12, 13]. The l.u.b. Lipschitz constant for a map f is de®ned as k f…u†ÿ f…v†k L‰ fŠˆ sup ; …19† kuÿ vk u6ˆv and the l.u.b. logarithmic norm (or l.u.b. Dahlquist constant) L‰I ‡ hfŠÿ 1 M‰ fŠˆ lim : …20† For a matrix A, we clearly have L‰AŠˆkAk and hence (20) and (12) coincide. With this extension one can study the growth of solutions and error growth of discretization methods for nonlinear ODEs y ˆ f…y†. Perhaps it is not so well known that these tools, in addition to the study of the growth of solutions and error growth on discretization methods, have also many other interesting applications. For example, given a matrix, as we have pointed out above, condition l‰AŠ > 0 is equivalent to the regularity of A. This ÿ1 is also true for any map, namely, if l‰ fŠ > 0, then f exists and ÿ1 ÿ1 L‰ f Šˆ l‰ fŠ (12). Moreover, as ÿl‰ fŠ M‰ fŠ L‰ fŠ a similar result can be given in terms of M‰ fŠ. More precisely, if ÿ1 ÿ1 ÿ1 M‰ fŠÿ< 0 then f exists and L‰ f ŠÿM‰ fŠ [12]. For unbounded operators the classical de®nitions cannot be used. However, in the case of Hilbert spaces, as we have an inner product norm we may use (14) and de®ne the logarithmic norm as, huÿ v; f…u†ÿ f…v†i M ‰ fŠˆ sup : …21† huÿ v; uÿ vi u6ˆv H STIFF INDEX-3 DAEs 247 In expression (21) it is not required f to be a bounded operator. For example, we may consider the 1D reaction±diffusion model equation [14] u ˆ u ‡ g…u† t xx in H ‰0; 1Š with boundary data u…t; 0†ˆ u…t; 1†ˆ 0. If hu; viˆ uv dx, 0 0 integration by parts yields hu; u iˆÿhu ; u iÿ hu; ui, where the last xx x x 2 2 2 inequality follows from Sobolev's lemma. Thus M ‰@ =@x Šˆÿ although the diffusion operator is unbounded. An extension of the concept of logarithmic norm for matrix pencils was done in [7]; the aim for this extension was to have a tool to study the qualitative behavior of the solution of linear variable coef®cient DAEs as well as their numerical approximations [4, 3]. The l.u.b. logarithmic norm for a matrix pencil is de®ned as L ‰A; A‡ BŠÿ 1 M ‰A; BŠˆ lim …22† !0 where V is a subspace such that V 6ˆ f0g and V \ Ker Aˆf0g ; …23† and kBxk L ‰A; BŠˆ max : …24† x2V ; x6ˆ0 kAxk When a subspace V satis®es (23), we say that it is an admissible subspace for L ‰A;Š. Observe that for V ˆ R , M ‰I;ÿBŠˆ M‰BŠ and L ‰I; BŠˆ L‰BŠ.In V V V this way, (24) and (22) are the extensions of (13) and (12) respectively. Finally, in Banach spaces, by the use of semi-inner product norms, logarithmic norms can also be de®ned. In this context, in [8] restricted Lipschitz constants and restricted logarithmic norms are introduced to study stability issues for problems of the form Ax ˆ f…x; t† with A singular. As the map f is not assumed to be bounded, index one partial differential algebraic equations are included in the study. After this introduction on l.u.b. / g.l.b. Lipschitz constants and l.u.b. / g.l.b. logarithmic norms constants if we go back to the DAE (1)±(2), we must realize that the ellipticity condition (3) is simply a restricted g.l.b logarithmic norm (18). Taking into account that the analysis of stability bounds to derive convergence for stiff ODEs is done in terms of M‰Š, and that m‰Š ˆ ÿM‰ÿŠ, 248 I. HIGUERAS we immediately may infer that the techniques used for ODEs can be transferred for DAEs. In some sense, it is simply a matter of language: whereas in the classical numerical analysis for ODEs the M‰Š language is more common, in the PDE ®eld, the m‰Š and l‰Š language is more usual. Once we have realized of this relationship, it will not be dif®cult to obtain stability and convergence results for (1)±(2). As our further aim is to study not only problems with this structure, we extend the g.l.b. Lipschitz constants l‰Š and g.l.b. logarithmic norms m‰Š to matrix pencils. To do so, we follow the lines of [7]. We should point out that in what concerns to m‰Š, more than a new extension, it is simply a translation of what is done in [7]. 3.1 G.l.b. Lipschitz Constants for Matrix Pencils Given a linear constant coef®cient DAEs with matrix pencil…A; B†, in order to get uniqueness for the solution of initial value problems, the pencil must be regular. This means that there exists  2 C such that det…†  A‡ B6ˆ 0. Otherwise, if det…†  A‡ B ˆ 0 for all , the pencil is called singular. Observe that if the pencil…A; B† is regular, then Ker A\ Ker Bˆf0g. Otherwise, for a v 2 Ker A\ Ker B,wehave … A‡ B†v ˆ 0 for any  , and the pencil is singular. For a regular pencil, det…†  A‡ B is a polynomial with degree less or equal to n. The complex values  such that det…†  A‡ Bˆ 0 are called ®nite eigenvalues of the pencil. Observe that not every pencil has eigenvalues. For 2 C, the vectors v 2 C such that … A‡ B†v ˆ 0 are called eigenvectors associated to the ®nite eigenvalue . Observe that for any eigenvector v associated to ®nite eigenvalues, we have Av 6ˆ 0. Given a matrix pencil…A; B†, we are looking for an extension of g.l.b. norm for a matrix, such that a similar expression to (3), whenever it exists, holds. In order to ®nd this value, we proceed in a similar way it is done for l.u.b. norms. If v is an eigenvector of the pencil associated to the eigenvalue ,we have, jjkAvkˆkBvk, and as Av 6ˆ 0, we can consider jjˆkBvk =kAvk.If V 6 Ker A, we can de®ne kBvk l ‰A; BŠˆ inf ; …25† v2VkAvk We give the following de®nition. De®nition 3.1 We say that a subspace V is an admissible subspace for l ‰A;Š if V 6 Ker A. V STIFF INDEX-3 DAEs 249 Observe that l ‰I; AŠˆ l‰AŠ: Observe too that if V 6ˆ f0g, then an IR admissible subspace for L ‰A;Š is also an admissible subspace for l ‰A;Š. V V The expressions l‰AŠ and L‰AŠ are related by (17), or in a equivalent way 0if R \ Ker A 6ˆ f0g l‰AŠˆ ÿ1 ÿ1 n L‰A Š if R \ Ker Aˆf0g It is straightforward to prove that a similar relationship exists between L ‰A; BŠ and l ‰A; BŠ. V V Proposition 3.1 For the g.l.b. Lipschitz constant, it holds 0if V \ Ker B 6ˆ f0g l ‰A; BŠˆ V ÿ1 L ‰B; AŠ if V \ Ker Bˆf0g Directly from the de®nition or using Proposition 3.1 and the properties of L ‰A; BŠ in [7], we can derive the following properties. Proposition 3.2 Properties of the g.l.b. Lipschitz constant (25): 1. kBvk l ‰A; BŠkAvk, 8v 2 V . j j 2. For any , 2 C, with 6ˆ 0, we have l ‰ A; BŠˆ l ‰A; BŠ. V V j j 3. l ‰A; B‡ CŠ l ‰A; BŠ‡ l ‰A; CŠ. V V V 4. If V \ Ker B 6ˆ f0g, then l ‰A; BŠˆ 0. Thus l ‰A; BŠˆ 0 does not imply V V B ˆ 0. The application P : L…R †! R ; B,!P …B†ˆ l ‰A; BŠ A;V A;V V is a seminorm. 5. If V contains an eigenvector associated with any eigenvalue  of …A; B†, then jj l ‰A; BŠ : 6. If V contains any eigenvector associated to any eigenvalue  of …A; B† such that jjˆ minfjj :  eigenvalue of …A; B†g, then i i jj l ‰A; BŠ : 7. Let P a non-singular matrix, then l ‰AP; BPŠˆ l ‰A; BŠ where V PV PV ˆfPx = x 2 Vg. In particular if A is non singular, then ÿ1 ÿ1 n n l ‰A; BŠˆ l ‰I; BA Šˆ l‰BA Š : IR IR 250 I. HIGUERAS 8. For any non-singular matrix P, in general l ‰PA; PBŠ6ˆ l ‰A; BŠ. V V 9. The de®nition depends on V . 10. If V  W , then l ‰A; BŠ l ‰A; BŠ. V W 11. l ‰A; BŠ L ‰A; BŠ. V V In the following table we remark the duality between L ‰A; BŠ and l ‰A; BŠ. V B V admissible: V \ Ker Aˆf0g V admissible: V 6 Ker A V  Ker B : L ‰A; BŠˆ 0 V \ Ker B 6ˆ f0g : l ‰A; BŠˆ 0 V V If we denote by  …A; B† the minimum eigenvalue of the pencil …A; B†,it min is easy to give a characterization of l ‰A; BŠ for the euclidean norm. Proposition 3.3 Consider the n s matrix S whose columns are a basis of V . Then p t t t t l ‰A; BŠˆ  …S A AS; S B BS† : V;2 min 3.2 G.l.b. Logarithmic Norms for Matrix Pencils With the de®nition of l ‰A; BŠ, we can now de®ne the g.l.b. Dahlquist constant of a matrix pencil as l ‰A; Aÿ BŠÿ 1 m ‰A; BŠˆ lim : …26† !0 Observe that for V ˆ R , m ‰I;ÿBŠˆ m‰BŠ. It can be proved that m ‰A; BŠˆÿM ‰A;ÿBŠ : V V From this relationship and using the properties of M ‰;Š in [7] or directly from the de®nition, we ®nd the following properties. Proposition 3.4 Properties of the g.l.b. logarithmic norm (26): 1. ÿl ‰A; BŠ m ‰A; BŠ l ‰A; BŠ. V V V 2. m ‰A; BŠ M ‰A; BŠ. V V 3. For any , in R, 6ˆ 0, j j m ‰ A; BŠˆ m ‰A; sgn …  † BŠ : V V j j STIFF INDEX-3 DAEs 251 4. m ‰A; B‡ CŠ m ‰A; BŠ‡ m ‰A; CŠ. V V V 5. If V contains an eigenvector associated with any eigenvalue  of …A; B†, then Re …† m ‰A; BŠ. 6. m ‰A; B‡ zAŠˆ m ‰A; BŠ‡ Re z for all z 2 C. V V 7. Ifkk is an inner product norm, then hAx;ÿBxi m ‰A; BŠˆ min : Ax6ˆ0 hAx; Axi x2V This means that l ‰A; BŠ is the greatest constant such that hAx;ÿBxi hAx; Axi for all x 2 V : 8. If A is invertible, and V ˆ R , then ÿ1 ÿ1 m ‰A; BŠˆ m ‰I ; BA Šˆ m‰ÿBA Š : V V n 9. kBxk max…m ‰A;ÿBŠ; m ‰A; BŠ†, for all x 2 V . V V 10. For regular matrices T and T, ÿ1 ÿ1 ~ ~ m …A; B†ˆ m …TA; TB†ˆ m …TAT ; TBT † V;T V TV and if V is T -invariant, ÿ1 ÿ1 ~ ~ m …A; B†ˆ m …TAT ; TBT † : V;T V Again for the euclidean norm there is a nice characterization. Proposition 3.5 Consider the n s matrix S whose columns are a basis of V . Then t t …B A‡ A B† t t t m ‰A; BŠˆ  S A AS;ÿS S : V;2 min 4 APPLICATION TO MECHANICAL SYSTEMS We consider the system 00 t u ‡ Au‡ B z ˆ f …27† Bu ˆ f …28† 2 252 I. HIGUERAS n m where u 2 R , z 2 R , A is an n n matrix such that m ‰AŠˆ > 0, i.e., Ker B hAu; ui kuk for all u 2 Ker B and B is an m n matrix such that l‰BŠˆ > 0, i.e., t m kB zk kzk for all u 2 R : As it has been pointed out in Remark 3.2, the condition imposed on B implies t t that l‰BB Š > 0 and hence the matrix BB is regular. Observe that we have dropped the condition on the continuity of A. 4.1 Perturbation Index We begin analyzing the system (28)±(28). Our aim is to obtain the perturbation index of the problem as well as to derive some information on the qualitative behavior of the solution. The DAE (27)±(28) is a particular case of the ones considered in [7]. The matrix pencil of the order 1 problem is 00 1 0 11 I 0 ÿI 0 ~ ~ @@ A @ AA A; B ˆ I ; A 0 B 0 B 00 The admissible subspace which contains the solution of the homogeneous problem is 2n‡m t ÿ1 V ˆf…u; v; z†2 R : u; v 2 Ker B ; zˆÿ…BB † BAug ; and the l.u.b. logarithmic norm was 0 ÿI ~ ~ M A; B ˆ M  I; V V A 0 0 I ˆ M ÿ A 0 ÿ1 2n t t with V ˆf…u; v†2 R : u; v 2 Ker Bg and :ˆ I ÿ B …BB † B. Observe that  is a projector onto Ker B, i.e.,  ˆ  and im  ˆ Ker B. Observe that with the euclidean inner product norm, the projector  is an orthogonal projector. STIFF INDEX-3 DAEs 253 ~ ~ A suf®cient condition for asymptotic stability of …u; v† is M ‰A; BŠ < 0, or in terms of the g.l.b. logarithmic norm, 0 ÿI m > 0 : A 0 If we have an inner product norm (Proposition 3.4-7)), 0 ÿI hAu; viÿhu; vi m ˆ min 2 2 u;v2Ker B A 0 kuk ‡kvk hAu; uiÿhu; ui min u2Ker B 2kuk 1 1 ˆ…† m ‰AŠÿ 1ˆ…† m ‰AŠÿ 1 ; Ker B Ker B 2 2 where we have used the orthogonality of . On the following we impose ~ ~ m A; B > 0 : …29† Observe that this condition implies that m ‰AŠ > 1 ; Ker B which is stronger than the Ker B-ellipticity of the map A (3). Recall that on the other hand we do not impose the continuity of A. Remark 4.1 In this part we deal with the ®nite dimensional index-3 DAE but having in mind PDAEs, i.e., _ the general operator index-3 DAE (7)±(8). Thus we must have care so that all the steps we do in the decoupling be also valid for that system. & First we are going to ®nd how the solution is decoupled in (27)±(28). It is well known that the differential index as well as the perturbation index of this DAE is 3. Recall that although in the PDAE framework the differential index has no sense, there is not any problem in de®ning the perturbation index. Thus we are going to decouple the DAE to prove that the perturbation index is 3 but using properties that are also valid for the operator DAE. This decoupling will be extremely useful for the analysis of the numerical solution. We consider consistent initial conditions, 0 0 u…t †ˆ u u…t †ˆ u 0 0 0 0 254 I. HIGUERAS satisfying 0 0 Bu ˆ f …t † ; Bu ˆ f …t † : 0 2 0 0 0 2 We begin with standard computations to obtain ÿ1 ÿ1 ÿ1 t t t 00 zˆÿ…BB † BAu‡…BB † Bf ÿ…BB † f …30† and substituting it into (27) we get ÿ1 00 t t 00 u ‡  Au ˆ f ‡ B …BB † f : …31† Any vector u can be written as u ˆ u‡…I ÿ †u. If we multiply (31) by I ÿ  and , and denote u :ˆ…I ÿ †u, u :ˆ u, we obtain 1 2 00 t t ÿ1 00 u ˆ B …BB † f 1 2 u ‡ Au ˆ f The initial condition is 0 0 Bu ˆ f …t † ; Bu ˆ f …t † ; 1;0 2 0 0 1;0 2 and thus the ®rst equation implies ÿ1 t t u ˆ B …BB † f : 1 2 From the second one we obtain ÿ1 00 t t u ‡ Au ˆ  f ÿ AB …BB † f : …32† 2 1 2 Recall that taking into account the origin of this ODE, when we deal with it, we must consider initial conditions satisfying …I ÿ †u …t †ˆ 0. 2 0 In order to study the qualitative the behavior of the solutions, it is suf®cient consider ODEs of the form (32). Thus we consider order 1 problem u ˆ v …33† v ˆÿAu‡  p …34† This is an ordinary differential equation with Im  as invariant subspace, i.e., if u…t †, v…t †2 Im , then u…t†, v…t†2 Im  for all t. To see this, we multiply 0 0 the second equation by I ÿ , to obtain …I ÿ †v ˆ 0, and hence …I ÿ † v…t†ˆ…I ÿ †v…t †. Hence if …I ÿ †v…t †ˆ 0, then …I ÿ †v…t†ˆ 0 and 0 0 therefore also …I ÿ †v …t†ˆ 0; equivalently, we obtain v…t†ˆ v…t† and STIFF INDEX-3 DAEs 255 u…t†ˆ u…t†. An important remark at this point is that the stability of (33)±(34) must be studied only on Im . This situation is not new, for example it also appears when the nonlinear index-1 DAE is decoupled [5]. Now the standard perturbations bounds for ODEs applied to (33)±(34) with invariant subspaces gives a…t† a…t† ÿa…s† k…u…t†; v…t††k  e k…u…t †; v…t †k‡ e e kp…s†k ds t  t 0 0 0 …35† ÿ1 t t where a…t†ˆ…tÿ t †m ‰A; BŠ, and p…t†ˆ f …t†ÿ  AB …BB † f …t†. 0 V 1 2 Therefore, for the solution of (27)±(28) we have kukku …t†k‡ku …t†k 2 1 a…t† 0 e…† ku…t †k‡ku…t †k 0 0 ÿ1 a…t† ÿa…s† t t ‡ e e kp…s†k ds‡kB …BB † f k ÿ1 ÿ1 ÿ1 t t t 00 kzkk…BB † BAkkuk‡k…BB † Bf ÿ…BB † f k In particular we observe that in the bounds we have f , and thus the perturbation index is 3. Condition (29) gives information on the asymp- totic stability of the problem. Recall that these bounds are also valid for (7)±(8). 4.2 Numerical Solutions by Runge±Kutta Methods In this section we show how, given (27)±(28), we can give convergence estimates independently of the stiffness of the matrix A, we will use the l.u.b. or g.l.b. constant instead. Recall that the problem we are considering is a simple lineal problem but we use techniques that can be used for for (7)±(8). We consider an s-stages Runge±Kutta method with coef®cients…b ;A†.We assume that the matrix A is regular and the method is stif¯y accurate, i.e., ÿ1 a ˆ b , i ˆ 1; .. . ; s. We denote t ˆ t ‡ c h, and ˆ…A † . As usual, si i ni nÿ1 i ij ij we assume that A l ˆ c and thus c ˆ 1. Given an approximation …u ; v ; z † of the solution at t , the new approximation …u ; v ; z † nÿ1 nÿ1 nÿ1 nÿ1 n n n at t ˆ t ‡ h is given by n nÿ1 u ˆ U ; v ˆ V ; z ˆ Z ; n ns n ns n ns 256 I. HIGUERAS where the internal stages U , V , Z , i ˆ 1; .. . ; s are obtained solving the ni ni ni system ‰UŠ ˆ V ni ni 0 t ‰VŠ ‡ AU ‡ B Z ˆ f ni ni 1i ni BU ˆ f i ˆ 1; .. . ; s ni 2i with the internal derivatives de®ned by s s X X ÿ ÿ 1 1 0 0 ‰UŠ ˆ U ÿ u ; ‰VŠ ˆ V ÿ v ij nj nÿ1 ij nj nÿ1 ni ni h h jˆ1 jˆ1 i ˆ 1; .. . ; s : We can state the following result. The proof is given in the Appendix. Theorem 4.1 Let an stif¯y accurate Runge±Kutta method …A; b † with positive weights be algebraically stable. Suppose that a positive diagonal matrix D exists such that …A† > 0. If q is the order stage of the method, then the method is convergent with order at least q for the variable u, q for the variable v, i.e., u , and qÿ 1 for the variable z . Thus for example, for the RadauIIA, we have that the orders are: s for the variable u, s for v, and sÿ 1 for z. The nonstiff order is 2sÿ 1, s and sÿ 1 respectively. Remark 4.2 Recall that Theorem 2.2, and therefore Theorem 4.1, is valid on the class F . Problems in this class must have smooth derivatives ([2, (7.2.2)]). Thus in particular it is not valid for highly oscillatory problems 00 2 like u ‡ w u ˆ 0. Roughly speaking, smoothness implies that the effective order stage is the one given by the algebraic conditions C…q†. & 5 CONCLUSIONS In this paper we have used logarithmic norms to obtain the perturbation index and convergence results for stif¯y accurate Runge±Kutta methods for a linear index-3 DAE coming from mechanical systems. This is the natural way if we want to make use of the experience of stiff ODEs [15, 2]. The decoupling for the linear index-3 DAE coming from mechanical systems is also valid for the general operator DAE (7)±(8), and thus we can de®ne the perturbation index STIFF INDEX-3 DAEs 257 also for these problems. Concerning to the numerical solution, the decoupling done is also valid for (7)±(8) but in order to ensure convergence, a detailed revision of the proofs of numerical methods for stiff ODEs [2] must be done. The open question is if properties on Kronecker products are also valid when a ®nite dimensional space and a Hilbert space are involved; these properties are implicitly used in the proofs. In any case, there is not any problem to deal with the ®nite dimensional case and, as it was expected, there is an order reduction. Taking into account that logarithmic norms are also de®ned for nonlinear maps and for unbounded operators, we conclude that this technique can also be used not only for nonlinear DAEs but also for DAEs expressed in terms of operators. This means that we have a tool to analyze solutions and numerical discretization in PDAEs, that together with the experience of theory of stiff ODEs [15, 2] may give good results. REFERENCES 1. Desoer, C. and Haneda, H.: The Measure of a Matrix as a Tool to Analyze Computer Algorithms for Circuit Analysis. IEEE Trans. Circuit Theory 19 (1972), pp. 480±486. 2. Dekker, K. and Verwer, J.G.: Stability of Runge±Kutta Methods for Stiff Nonlinear Differential Equations. North-Holland, Amsterdam, 1984. 3. Garcõ Âa-Celayeta, B.: Stability for Differential Algebraic Equations. PhD. Universidad Pu  blica de Navarra (1998). 4. Garcõ Âa-Celayeta, B. and Higueras, I.: Runge±Kutta Methods for DAEs. A New Approach. J. Comput. Appl. Math. 111 (1999), pp. 49±61. 5. Griepentrog, E. and Ma Èrz, R.: Differential Algebraic Equations and Their Numerical Treatment. Teubner Texte zur Mathematik 88, Leipzig, 1986. 6. Hairer, L., Lubich, Ch. and Roche, M.: The Numerical Solution of Differential Algebraic Equations by Runge±Kutta Methods. Lectures Notes in Mathematics 1409, Springer, Berlin, 1989. 7. Higueras, I. and Garcõa-Celayeta, B.: Logarithmic Norms for Matrix Pencils. SIAM J. Matrix Anal. Appl., 20 (1999), pp. 646±666. 8. Higueras, I. and So È derlind, G.: Logarithmic Norms and Nonlinear DAE Stability. Preprint Opto Matema Âtica e Informa Âtica. Universidad Pu  blica de Navarra Nr. 5±2001, 2001. 9. Jay, L.O.: Me Âthodes du type Runge±Kutta pour des e Âquations diffe Ârentielles alge Âbriques  Á d'index 3 avec des applications aux systemes Hamiltoniens. PhD. University of Geneve (1994). 10. Lubich, Ch.: Integration of Stiff Mechanical Systems by Runge±Kutta Methods. ZAMP,44 (1993), pp. 1022±1053. 11. Simeon, B.: Numerische Simulation gekoppelter Systeme von partiellen und differential- algebraischen Gleichungen in der Mehrk operdynamic. Habilitationsschrift, Universitat Karlsruhe, 1999. 258 I. HIGUERAS 12. So È derlind, G.: On Nonlinear Difference and Differential Equations. BIT 24 (1984), pp. 667±680. 13. So È derlind, G.: Bounds on Nonlinear Operators in Finite-dimensional Banach Spaces. Numer. Math. 50 (1986), pp. 27±44. 14. Soderlind, G.: Seminar Course Stability Concepts in Numerical Analysis. University of Lund, 1999. 15. Strom, T.: On Logarithmic Norms. SIAM J. Numer. Anal. 12 (1975), pp. 741±753. APPENDIX Proof of Theorem 4.1. In vectorial form the method can be written as ‰UŠ ˆ V …36† ‰VŠ ‡D U ‡D t Z ˆ f …37† A n B n 1 D U ˆ f i ˆ 1; .. . ; s …38† B n 2 where we have denoted U ˆ…U ; .. . ; U †, and in a similar way V , Z , n n1 ns n n 0 0 0 0 0 ‰UŠ ˆ…‰UŠ ; .. . ;‰UŠ † and in a similar way‰VŠ ,‰ZŠ ,D ˆ I A, with A s n n1 ns n n the Kronecker product [2]. In vectorial form, the internal derivatives can be expressed as ÿ1 ‰UŠ ˆ …A I†…† U ÿ l u ; n nÿ1 0 ÿ1 ‰VŠ ˆ …A I†…† V ÿ l v ; n nÿ1 with lˆ…1; .. . ; 1†2 R . In terms of the internal stages, (36)±(38) is ÿ1 …A I†…† U ÿ l u ˆ hV n nÿ1 n ÿ1 …A I†…† V ÿ l v ‡ hD U ‡ hD t Z ˆ hf n nÿ1 A n B n 1 D U ˆ f : B n 2 As usual, in order to study convergence, we consider the perturbed system ÿ1 …A I†…† U ÿ l u ˆ hV ‡ R n nÿ1 n 1 ÿ1 …A I†…† V ÿ l v ‡ hD U ‡ hD t Z ˆ hf ‡ R n nÿ1 A n B n 1 2 D U ˆ f ‡ R B n 2 3 STIFF INDEX-3 DAEs 259 We denote E ˆ U ÿ U , and in a similar way E , E , and e ˆ u ÿu , u;n n n v;n z;n u;n n n and in a similar way e , e , to obtain v;n z;n ÿ1 …A I† E ÿ l e ˆ hE ; n‡ R …39† u;n u;nÿ1 v 1 ÿ1 …A I† E ÿ l e ‡ hD E ‡ hD E ˆ R …40† v;n v;nÿ1 A u;n B z;n 2 D E ˆ R …41† B u;n 3 We multiply (40) by D , to obtain ÿ1 …A I†D E ÿD …l e † ‡ hD E ‡ hD E ˆD R B v;n B v;nÿ1 BA u;n BB z;n B 2 As BB is regular, ÿ1 hE ˆÿ…A I†D ÿ1 E ÿD ÿ1 …l e † z;n t v;n t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R t t u;n 2 …BB † BA …BB † B Multiplying (39) by D , we obtain ÿ1 …A I†D E ÿD …l e † ˆ hD E ‡D R B u;n B u;nÿ1 B v;n B 1 and hence 1 1 ÿ1 D E ˆ …A I†D E ÿD …l e † ÿ D R B v;n B u;n B u;nÿ1 B 1 h h 1 1 ÿ1 ˆ …A I† R ÿD …l e † ÿ D R 3 B u;nÿ1 B 1 h h where we have used (41). Thus, ÿ2 hE ˆÿ …A I†D ÿ1 R z;n t 3 …BB † ÿ2 ‡ …A I†D ÿ1 …l e † u;nÿ1 …BB † B ÿ1 ÿ1 ÿ1 ÿ1 ‡ …A I†D R ‡…A I†D …l e † t 1 t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R …42† t u;n t 2 …BB † BA …BB † B Finally, substituting it into (40), we get ÿ1 ÿ1 …A I†E ÿ…A I†D …l e †‡ hD E v;n  v;nÿ1 A u;n ÿ2 ÿ1 ˆD R ‡ …A I†D R 2 t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h 260 I. HIGUERAS Together with (39), we have a system of the form ÿ1 …A I†E ÿ hE ˆ T …43† u;n v;n 1;n ÿ1 hD E ‡…A I†E ˆ T …44† A u;n v;n 2;n where ÿ1 T ˆ…A I†…l e †‡ R 1;n u;nÿ1 1 ÿ1 ÿ2 T ˆ…A I†D …l e †‡D R ‡ …A I†D ÿ1 R 2;n  v;nÿ1  2 t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h Multiplying (43)±(44) by  and I ÿ , we obtain ÿ1 …A I†D E ÿ hD E ˆD T …45† u;n  v;n  1;n ÿ1 hD E ‡…A I†D E ˆD T …46† A u;n  v;n  2;n and ÿ1 …A I†D E ÿ hD E ˆD T …Iÿ† u;n …Iÿ† v;n …Iÿ† 1;n ÿ1 …A I†D E ˆD T …Iÿ† v;n …Iÿ† 2;n respectively, with ÿ1 D T ˆ…A I†D …l e †‡D R 1;n  u;nÿ1  1 ÿ1 D T ˆ…A I†D …l e †‡D R 2;n  v;nÿ1  2 and ÿ1 D T ˆ…A I†D …l e †‡D R …Iÿ† 1;n …Iÿ† u;nÿ1 …Iÿ† 1 ÿ2 D T ˆ …A I†D ÿ1 R …Iÿ† 2;n t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h Thus, for the I ÿ  part of the error, we have D E ˆ h…A I†D T ‡…A I†D T u;n 2;n 1;n …Iÿ† …Iÿ† …Iÿ† D E ˆ…A I†D T …Iÿ† v;n …Iÿ† 2;n STIFF INDEX-3 DAEs 261 and substituting, D T and D T we obtain …Iÿ† 1;n …Iÿ† 2;n D E ˆD ÿ1 R …47† t t …Iÿ† u;n 3 B …BB † ÿ1 D E ˆ …A I†D ÿ1 R …48† v;n t t 3 …Iÿ† B …BB † 1 1 ÿ1 ÿ …A I†D …l e †ÿ D R …49† …Iÿ† u;nÿ1 …Iÿ† 1 h h In particular, if R ˆ 0, thenD E ˆ 0, and thus…I ÿ †e ˆ 0. In this 3 u;n u; n …Iÿ† case we obtain D E ˆ 0 …50† …Iÿ† u;n D E ˆÿ D R …51† …Iÿ† v;n …Iÿ† 1 In any case, the errors in these components are not propagated in the next steps. …2† Denoting E :ˆD E for the  part of the error, we have a system of the ;n ;n form ÿ1 …2† …2† …A I†E ÿ hE ˆ  …52† 1;n u;n v;n ÿ1 …2† …2† hD E ‡…A I†E ˆ  …53† A 2;n u;n v;n with ˆD T ;  ˆD T ÿ hD E : 1;n  1;n 2;n  2;n A…Iÿ† u;n Observe that this is the numerical solution with the Runge±Kutta method of the ODE, u 0 ÿI u 0 2 2 ‡ ˆ v A 0 v 0 2 2 with Im  ˆ Ker B as invariant subspace and (29) 0 I 0 ÿI ÿM ˆ m > 0 : V V ÿ A 0  A 0 ÿ1 k ‡1 t t k ‡1 1 2 If R ˆ #…h †, and R ÿ A…I ÿ †B …BB † R ˆ #…h †, the order 1 2 3 stage we have to consider in (52)±(53) is k ˆ minfk ; k g. Thus we can apply 1 2 Theorem 2.2 of the stiff theory of ODEs to obtain B-convergence of order at least the order stage, i.e., k. 262 I. HIGUERAS q‡1 q‡1 If the RK method has order stage q, then R ˆ #…h †, R ˆ #…h †. 1 2 For the exact numerical solution, R ˆ 0. Thus for the  part of the solution we have order at least q. As for the I ÿ  part of the solution we had (50)±(51), D E ˆ 0 …Iÿ† u;n D E ˆÿ D R ; …Iÿ† v;n …Iÿ† 1 and for the error in the component z, we had (42), ÿ1 ÿ1 hE ˆ …A I†D ÿ1 R ‡…A I†D ÿ1 …l e † z;n t 1 t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R t u;n t 2 …BB † BA …BB † B http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

Numerical Methods for Stiff Index-3 DAEs

Download PDF
24 pages

Loading next page...
 
/lp/taylor-francis/numerical-methods-for-stiff-index-3-daes-zYGNHTyYZv

References (17)

Publisher
Taylor & Francis
Copyright
Copyright Taylor & Francis Group, LLC
ISSN
1744-5051
eISSN
1387-3954
DOI
10.1076/mcmd.7.2.239.3648
Publisher site
See Article on Publisher Site

Abstract

Mathematical and Computer Modelling of Dynamical Systems 1387-3954/01/0702-239$16.00 2001, Vol. 7, No. 2, pp. 239±262 Swets & Zeitlinger I. HIGUERAS ABSTRACT Space semidiscretization of PDAEs, i.e. coupled systems of PDEs and algebraic equations, give raise to stiff DAEs and thus the standard theory of numerical methods for DAEs is not valid. As the study of numerical methods for stiff ODEs is done in terms of logarithmic norms, it seems natural to use also logarithmic norms for stiff DAEs. In this paper we show how the standard conditions imposed on the PDAE and the semidiscretized problem are formally the same if they are expressed in terms of logarithmic norms. To study the mathematical problem and their num- erical approximations, this link between the standard conditions and logarithmic norms allow us to use for stiff DAEs techniques similar to the ones used for stiff ODEs. The analysis is done for problems which appear in the context of elastic multibody systems, but once the tools, i.e., logarithmic norms, are developed, they can also be used for the analysis of other PDAEs / DAEs. Keywords: DAEs, elastic multibody systems, logarithmic norms, Runge±Kutta method, stiff. 1 INTRODUCTION Nowadays there is an increasing interest in the analysis and ef®cient numerical solution of coupled systems of partial differential equations (PDEs) and algebraic equations. This kind of problems appear for example in models of constrained mechanical systems including rigid and elastic bodies [11]. In some ®elds, these coupled systems are known as partial differential-algebraic systems (PDAEs). The aim of this paper is to analyze and obtain valid bounds for the numer- ical solution of PDAEs / DAEs which appear in the context of mechanical systems. We will deal with problems of the form 00 0 u ‡ Au‡ B  ˆ l …1† Bu ˆ m …2† Departamento de Matema Âtica e Informa Âtica, Universidad Pu  blica de Navarra, 31006 Pamplona, Spain. E-mail: higueras@unavarra.es 240 I. HIGUERAS 0 0 with A and B time independent linear maps, A : V ! V , B : V ! Q , V and Q 0 0 0 are Hilbert spaces, V and Q are their dual spaces respectively, and B denotes the dual map of B. The map A is such that there exists a constant > 0 with hAu; ui ;kuk for all u 2 Ker B ; …3† and B such that there exists a constant > 0 with kB k kk for all  2 Q : …4† n m 0 t For the ®nite dimensional case, we will take V ˆ R , Q ˆ R , and B ˆ B . Recall that we have in mind not only DAEs in the standard case but also DAEs in the general operator sense. The reason to consider only these simple systems is that they cover both PDAEs formulated as operator DAEs and DAEs obtained from a PDAE by a semidiscretization process in space. As these problems are stiff or may contain unbounded operators, the standard theory of numerical methods for DAEs [6, 9] is not valid and a new analysis must be done. Taking into account that in the theory of numerical methods for stiff ODEs the tool used is the logarithmic norm of a map (more precisely, the l.u.b. logarithmic norm), it seems natural to use logarithmic norms also for stiff DAEs. This paper shows how the standard conditions imposed on the operator formulation of the PDAE and the semidiscretized problem are formally the same if they are expressed in terms of logarithmic norms. To study the mathematical problem and the numerical approximations, this link between the standard conditions and logarithmic norms allows us to use for stiff DAEs similar techniques to the ones used for stiff ODEs. In Section 1, we give the PDAE formulation as operational DAEs and the semidiscretization process, and show that (1)±(2) together with the conditions (3)±(4) cover both problems. The section ends with a brief outline of numerical methods for stiff ODEs. In Section 2, devoted to logarithmic norms, we summarize the main de®nitions, extensions and results on this topic. The origin for this paper are problems which appear in the context of elastic multibody systems (1)±(2), but we also have in mind more general DAEs. This is the reason to extend the concept of g.l.b. logarithmic norm of a matrix to a matrix pencil following the lines of the one done in [7] for the l.u.b. logarithmic norm of a matrix. Once we have constructed the tools, we can use them to analyze the system (1)±(2) and its numerical solution with stif¯y accurate Runge±Kutta methods. This is done in Section 4. Finally, in Section 5 some conclusions and forthcoming work is pointed out. STIFF INDEX-3 DAEs 241 2 MATHEMATICAL MODELS AND NUMERICAL METHODS BACKGROUND Mathematical modeling may lead to coupled systems of partial differential equations and algebraic equations. In [11], in the context of elastic multibody systems, the following PDAE in weak formulation is studied, ``Find …u…; t†;…; t†† 2 V  Q such that hu ; vi‡ a…u; v†‡ b…v;†ˆhl; vi for all v 2 V …5† b…u;†ˆhm;i for all  2 Q '' ; …6† where V and Q are two Hilbert spaces, and a : V  V ! R, and b : V  Q ! R are two continuous bilinear forms. If we denote V ˆfv 2 V : b…v;†ˆ 0 for all  2 Qg ; it is assumed that the bilinear form a…;† is V -elliptic, i.e., there exists a constant > 0 such that a…u; u† kuk for all u 2 V : Moreover, the bilinear form b…;† is assumed continuous, i.e., there exists a constant  such that ka…u; v†k  kukkvk for all u 2 V : A standard condition for the bilinear form b…;† is the LBB-condition, i.e., there exists a constant > 0 such that b…v;† inf sup  : 2Q kvk kk v2V V Q 0 0 As it is pointed out in [11], we can consider the operators A : V ! V , with V the dual space, de®ned as hAu; viˆ a…u; v†, B : V ! Q , de®ned as 0 0 hBv;iˆ b…v;† and its dual operator B : Q ! V , de®ned as hB ; viˆ b…v;†. In this way, (5)±(6) is 00 0 u ‡ Au‡ B  ˆ l …7† Bu ˆ m …8† In terms of the map B, we have that V ˆ Ker B. The V -ellipticity of the 0 0 bilinear form a…;† can be written as hAu; ui kuk for all u 2 Ker B 242 I. HIGUERAS with > 0, the continuity of the bilinear form a…;† as kAuk kuk for all u 2 V ; and the LBB-condition as kB k kk for all  2 Q with > 0. Following the PDE approach, a method to solve numerically PDAEs consists in a semidiscretization in space followed by a time integration. In order to semidiscretize in space the equations (5)±(6), ®nite-dimensional subspaces V  V and Q  Q are considered, and solutions u …; t†2 V and  …; t†2 Q of the form q n X X u …x; t†ˆ v …x†q…t†; …x; t†ˆ …x†…t† i i  i i iˆ1 iˆ1 are constructed. If given a PDE the semidiscretization process gives an ODE system, for PDAEs, the problem obtained is a DAE system. The time dependent functions q…t†, …t† must solve the DAE i j 00 t M q …t†‡ A q…t†‡ B …t†ˆ l …t†…9† B q…t†ˆ m …t†…10† The mass matrix M and the matrix A are symmetric positive and semide®nite positive respectively. In (9)±(10) the matrices are de®ned in terms of the inner products in the Hilbert space and the bilinear forms, …M † ˆhv ; vi…A † ˆ a…v ; v†…B † ˆ b…v ;†: i j  i j  j i ij ij ij A compatibility condition is imposed on the discretization, 2 Q with b…v;†ˆ 0 for all v 2 V )  ˆ 0 which in terms of B is kB k kk with > 0 : With this condition the matrix B has full row rank and therefore the linear DAE (9)±(10) has index 3. To analyze convergence, the following subspace is considered W ˆf v 2 V : b…v;†ˆ 0 for all  2 Q g ; STIFF INDEX-3 DAEs 243 and it is imposed that the bilinear form a…;† is W -elliptic. Under these conditions, the matrix A satis®es hA q; qi kqk for all q 2 Ker B with > 0. The next step is to integrate (9)±(10) with some DAE method. We can expect that the nature of these DAEs be similar to the ODEs which result form the PDEs case: stiff problems. Therefore, the standard bounds which ensure certain order of convergence for index-3 DAEs [6, 9] are not valid anymore. Some authors [10] have studied the numerical solution of stiff mechanical systems using some model problems in which the stiffness parameter appears explicitly. As it was expected from the knowledge of stiff ODEs, depending on the stiffness parameter, an order reduction may appear. In the theory of numerical methods for stiff ODE problems, in order to analyze the effective order, the concept of B-convergence is introduced [2]. Basically the order of B-convergence is the order of the method independently of the stiffness of the problem. Other concepts used in this framework are the B-consistency and C-stability, which again are consistency and stability independently of the stiffness. Now the result is that B-consistency and C- stability imply B-convergence [2, Th. 7.2.5], and conditions to obtain C- stability and B- consistency must be studied. The class of ODEs for which this theory is developed is y ˆ f…t; y† such that h f…t; y †ÿ f…t; y †; y ÿ y i ky ÿ y k ; …11† 1 2 1 2 1 2 where h;i is an inner product norm in R . The class F is the class of problems satisfying (11) and such that their solutions are smooth, in the sense that the derivatives are bounded by constants independent of any quantity which may be in¯uenced by the stiffness [2, (7.2.2)]. The constant  is also known in the literature as the logarithmic norm of f , or the l.u.b. Dalhquist constant [2, 13]. When implicit Runge±Kutta methods are used, some nonlinear systems have to be used. In order to obtain existence and uniqueness of solution for these systems, the case < 0, is extremely favorable. Theorem 2.1 [2, Cor. 5.3.13] If a Runge±Kutta method is algebraically stable and irreducible, then the system of algebraic equations has unique solution for any problem from F with < 0. The important point in this result is that there is not any stepsize restriction. If > 0, a similar result is obtained with a stepsize restriction. 244 I. HIGUERAS Given a Runge±Kutta method with matrix coef®cients A, and in [2, (5.1.21)] the function hA ; i …A† ˆ min 6ˆ0 h;i is de®ned, where h;i denotes the inner product ~ ~ ~ h; i :ˆ d   ;  2 R ; d > 0 ; j ˆ 1; .. . ; s : j j j j jˆ1 In terms of this function, the following result is proved. The order of convergence is obtained from Theorem 7.3.3 in the same reference. Theorem 2.1 [2, Th. 7.4.4] Let a Runge±Kutta method…A; b † with positive weights be algebraically stable. Suppose that a positive diagonal matrix D exists such that …A† > 0. If q is the order stage of the method, then the method is B-convergent with order at least q on F . Again, the case < 0 is favorable in the sense that all the stability bounds involved to prove this result hold without any stepsize restriction. As a corollary, we obtain that the Gauss (order stage s), RadauIA (order stage sÿ 1), RadauIIA (order stage s), and the two stage Lobatto IIIC (order stage sÿ 1) methods are B-convergent. B-convergence of some other diagonally implicit and singly implicit methods are also proved in [2, Cor. 7.4.6, 7.4.7]. The numerical experiments shows that the lower bound for the order of convergence, i.e., the order stage of the method, can be reached. Recall that these proofs are done in terms of logarithmic norms and inner product norms, and they only use the standard properties of these concepts. 3 LOGARITHMIC NORMS The concept of logarithmic norm for a matrix was introduced in 1958 by Dahlquist and Lozinskij as a tool to study the growth of solutions to ODEs and the error growth in discretization methods for their approximate solution. For details see [2]. For a matrix A, the logarithmic norm (more precisely l.u.b. logarithmic norm) is de®ned by L‰I ‡  AŠÿ 1 M‰AŠˆ lim : …12† !0 STIFF INDEX-3 DAEs 245 The norm here is generally assumed to be an l.u.b. norm (or an l.u.b. Lipschitz constant), kAxk L‰AŠˆ max : …13† x6ˆ0 kxk Properties of norms and logarithmic norms can be found in [2, 1]. In particular, if the used norm is an inner product norm, we have the equivalent de®nition hAu; ui M‰AŠˆ max : …14† u6ˆ0 hu; ui In a similar way we can de®ne the g.l.b. logarithmic norms as l‰I ‡  AŠÿ 1 m‰AŠˆ lim ; …15† !0 where l‰Š is the g.l.b. Lipschitz constant kAxk l‰AŠˆ min : …16† x6ˆ0 kxk It can be proved that for the euclidean norm, p l ‰AŠˆ minf  :  eigenvalue of A Ag 2 i i ˆ minf :  singular value of Ag : i i Observe that from the de®nition l‰AŠ minfjj : det…I ÿ A†ˆ 0g : Moreover it can be proved that 0 if A singular l‰AŠˆ …17† ÿ1 ÿ1 L‰A Š if A regular and that m‰AŠˆÿM‰ÿAŠ. Thus if the norm is an inner product norm, hAu; ui m‰AŠˆ min : …18† u6ˆ0 hu; ui With the l.u.b. and g.l.b. logarithmic norms of a matrix, given the ODE y ˆ A…t†y, one can derive bounds for the solution and therefore give suf®cient 246 I. HIGUERAS conditions for stability. Some well-known theorems in stability theory of ODEs are easily proven with the help of logarithmic norms [15]. Remark 3.1 In the literature, usually these concepts are de®ned for square matrices but they can also be de®ned for rectangular matrices. & Remark 3.2 Observe that condition (4) can be expressed in terms of the g.l.b. Lipschitz constant l‰BŠˆ > 0. Hence we have t 2 l‰BBŠˆ > 0 ; and thus BB is a regular matrix. & Logarithmic norms for nonlinear maps are considered in [2, 12, 13]. The l.u.b. Lipschitz constant for a map f is de®ned as k f…u†ÿ f…v†k L‰ fŠˆ sup ; …19† kuÿ vk u6ˆv and the l.u.b. logarithmic norm (or l.u.b. Dahlquist constant) L‰I ‡ hfŠÿ 1 M‰ fŠˆ lim : …20† For a matrix A, we clearly have L‰AŠˆkAk and hence (20) and (12) coincide. With this extension one can study the growth of solutions and error growth of discretization methods for nonlinear ODEs y ˆ f…y†. Perhaps it is not so well known that these tools, in addition to the study of the growth of solutions and error growth on discretization methods, have also many other interesting applications. For example, given a matrix, as we have pointed out above, condition l‰AŠ > 0 is equivalent to the regularity of A. This ÿ1 is also true for any map, namely, if l‰ fŠ > 0, then f exists and ÿ1 ÿ1 L‰ f Šˆ l‰ fŠ (12). Moreover, as ÿl‰ fŠ M‰ fŠ L‰ fŠ a similar result can be given in terms of M‰ fŠ. More precisely, if ÿ1 ÿ1 ÿ1 M‰ fŠÿ< 0 then f exists and L‰ f ŠÿM‰ fŠ [12]. For unbounded operators the classical de®nitions cannot be used. However, in the case of Hilbert spaces, as we have an inner product norm we may use (14) and de®ne the logarithmic norm as, huÿ v; f…u†ÿ f…v†i M ‰ fŠˆ sup : …21† huÿ v; uÿ vi u6ˆv H STIFF INDEX-3 DAEs 247 In expression (21) it is not required f to be a bounded operator. For example, we may consider the 1D reaction±diffusion model equation [14] u ˆ u ‡ g…u† t xx in H ‰0; 1Š with boundary data u…t; 0†ˆ u…t; 1†ˆ 0. If hu; viˆ uv dx, 0 0 integration by parts yields hu; u iˆÿhu ; u iÿ hu; ui, where the last xx x x 2 2 2 inequality follows from Sobolev's lemma. Thus M ‰@ =@x Šˆÿ although the diffusion operator is unbounded. An extension of the concept of logarithmic norm for matrix pencils was done in [7]; the aim for this extension was to have a tool to study the qualitative behavior of the solution of linear variable coef®cient DAEs as well as their numerical approximations [4, 3]. The l.u.b. logarithmic norm for a matrix pencil is de®ned as L ‰A; A‡ BŠÿ 1 M ‰A; BŠˆ lim …22† !0 where V is a subspace such that V 6ˆ f0g and V \ Ker Aˆf0g ; …23† and kBxk L ‰A; BŠˆ max : …24† x2V ; x6ˆ0 kAxk When a subspace V satis®es (23), we say that it is an admissible subspace for L ‰A;Š. Observe that for V ˆ R , M ‰I;ÿBŠˆ M‰BŠ and L ‰I; BŠˆ L‰BŠ.In V V V this way, (24) and (22) are the extensions of (13) and (12) respectively. Finally, in Banach spaces, by the use of semi-inner product norms, logarithmic norms can also be de®ned. In this context, in [8] restricted Lipschitz constants and restricted logarithmic norms are introduced to study stability issues for problems of the form Ax ˆ f…x; t† with A singular. As the map f is not assumed to be bounded, index one partial differential algebraic equations are included in the study. After this introduction on l.u.b. / g.l.b. Lipschitz constants and l.u.b. / g.l.b. logarithmic norms constants if we go back to the DAE (1)±(2), we must realize that the ellipticity condition (3) is simply a restricted g.l.b logarithmic norm (18). Taking into account that the analysis of stability bounds to derive convergence for stiff ODEs is done in terms of M‰Š, and that m‰Š ˆ ÿM‰ÿŠ, 248 I. HIGUERAS we immediately may infer that the techniques used for ODEs can be transferred for DAEs. In some sense, it is simply a matter of language: whereas in the classical numerical analysis for ODEs the M‰Š language is more common, in the PDE ®eld, the m‰Š and l‰Š language is more usual. Once we have realized of this relationship, it will not be dif®cult to obtain stability and convergence results for (1)±(2). As our further aim is to study not only problems with this structure, we extend the g.l.b. Lipschitz constants l‰Š and g.l.b. logarithmic norms m‰Š to matrix pencils. To do so, we follow the lines of [7]. We should point out that in what concerns to m‰Š, more than a new extension, it is simply a translation of what is done in [7]. 3.1 G.l.b. Lipschitz Constants for Matrix Pencils Given a linear constant coef®cient DAEs with matrix pencil…A; B†, in order to get uniqueness for the solution of initial value problems, the pencil must be regular. This means that there exists  2 C such that det…†  A‡ B6ˆ 0. Otherwise, if det…†  A‡ B ˆ 0 for all , the pencil is called singular. Observe that if the pencil…A; B† is regular, then Ker A\ Ker Bˆf0g. Otherwise, for a v 2 Ker A\ Ker B,wehave … A‡ B†v ˆ 0 for any  , and the pencil is singular. For a regular pencil, det…†  A‡ B is a polynomial with degree less or equal to n. The complex values  such that det…†  A‡ Bˆ 0 are called ®nite eigenvalues of the pencil. Observe that not every pencil has eigenvalues. For 2 C, the vectors v 2 C such that … A‡ B†v ˆ 0 are called eigenvectors associated to the ®nite eigenvalue . Observe that for any eigenvector v associated to ®nite eigenvalues, we have Av 6ˆ 0. Given a matrix pencil…A; B†, we are looking for an extension of g.l.b. norm for a matrix, such that a similar expression to (3), whenever it exists, holds. In order to ®nd this value, we proceed in a similar way it is done for l.u.b. norms. If v is an eigenvector of the pencil associated to the eigenvalue ,we have, jjkAvkˆkBvk, and as Av 6ˆ 0, we can consider jjˆkBvk =kAvk.If V 6 Ker A, we can de®ne kBvk l ‰A; BŠˆ inf ; …25† v2VkAvk We give the following de®nition. De®nition 3.1 We say that a subspace V is an admissible subspace for l ‰A;Š if V 6 Ker A. V STIFF INDEX-3 DAEs 249 Observe that l ‰I; AŠˆ l‰AŠ: Observe too that if V 6ˆ f0g, then an IR admissible subspace for L ‰A;Š is also an admissible subspace for l ‰A;Š. V V The expressions l‰AŠ and L‰AŠ are related by (17), or in a equivalent way 0if R \ Ker A 6ˆ f0g l‰AŠˆ ÿ1 ÿ1 n L‰A Š if R \ Ker Aˆf0g It is straightforward to prove that a similar relationship exists between L ‰A; BŠ and l ‰A; BŠ. V V Proposition 3.1 For the g.l.b. Lipschitz constant, it holds 0if V \ Ker B 6ˆ f0g l ‰A; BŠˆ V ÿ1 L ‰B; AŠ if V \ Ker Bˆf0g Directly from the de®nition or using Proposition 3.1 and the properties of L ‰A; BŠ in [7], we can derive the following properties. Proposition 3.2 Properties of the g.l.b. Lipschitz constant (25): 1. kBvk l ‰A; BŠkAvk, 8v 2 V . j j 2. For any , 2 C, with 6ˆ 0, we have l ‰ A; BŠˆ l ‰A; BŠ. V V j j 3. l ‰A; B‡ CŠ l ‰A; BŠ‡ l ‰A; CŠ. V V V 4. If V \ Ker B 6ˆ f0g, then l ‰A; BŠˆ 0. Thus l ‰A; BŠˆ 0 does not imply V V B ˆ 0. The application P : L…R †! R ; B,!P …B†ˆ l ‰A; BŠ A;V A;V V is a seminorm. 5. If V contains an eigenvector associated with any eigenvalue  of …A; B†, then jj l ‰A; BŠ : 6. If V contains any eigenvector associated to any eigenvalue  of …A; B† such that jjˆ minfjj :  eigenvalue of …A; B†g, then i i jj l ‰A; BŠ : 7. Let P a non-singular matrix, then l ‰AP; BPŠˆ l ‰A; BŠ where V PV PV ˆfPx = x 2 Vg. In particular if A is non singular, then ÿ1 ÿ1 n n l ‰A; BŠˆ l ‰I; BA Šˆ l‰BA Š : IR IR 250 I. HIGUERAS 8. For any non-singular matrix P, in general l ‰PA; PBŠ6ˆ l ‰A; BŠ. V V 9. The de®nition depends on V . 10. If V  W , then l ‰A; BŠ l ‰A; BŠ. V W 11. l ‰A; BŠ L ‰A; BŠ. V V In the following table we remark the duality between L ‰A; BŠ and l ‰A; BŠ. V B V admissible: V \ Ker Aˆf0g V admissible: V 6 Ker A V  Ker B : L ‰A; BŠˆ 0 V \ Ker B 6ˆ f0g : l ‰A; BŠˆ 0 V V If we denote by  …A; B† the minimum eigenvalue of the pencil …A; B†,it min is easy to give a characterization of l ‰A; BŠ for the euclidean norm. Proposition 3.3 Consider the n s matrix S whose columns are a basis of V . Then p t t t t l ‰A; BŠˆ  …S A AS; S B BS† : V;2 min 3.2 G.l.b. Logarithmic Norms for Matrix Pencils With the de®nition of l ‰A; BŠ, we can now de®ne the g.l.b. Dahlquist constant of a matrix pencil as l ‰A; Aÿ BŠÿ 1 m ‰A; BŠˆ lim : …26† !0 Observe that for V ˆ R , m ‰I;ÿBŠˆ m‰BŠ. It can be proved that m ‰A; BŠˆÿM ‰A;ÿBŠ : V V From this relationship and using the properties of M ‰;Š in [7] or directly from the de®nition, we ®nd the following properties. Proposition 3.4 Properties of the g.l.b. logarithmic norm (26): 1. ÿl ‰A; BŠ m ‰A; BŠ l ‰A; BŠ. V V V 2. m ‰A; BŠ M ‰A; BŠ. V V 3. For any , in R, 6ˆ 0, j j m ‰ A; BŠˆ m ‰A; sgn …  † BŠ : V V j j STIFF INDEX-3 DAEs 251 4. m ‰A; B‡ CŠ m ‰A; BŠ‡ m ‰A; CŠ. V V V 5. If V contains an eigenvector associated with any eigenvalue  of …A; B†, then Re …† m ‰A; BŠ. 6. m ‰A; B‡ zAŠˆ m ‰A; BŠ‡ Re z for all z 2 C. V V 7. Ifkk is an inner product norm, then hAx;ÿBxi m ‰A; BŠˆ min : Ax6ˆ0 hAx; Axi x2V This means that l ‰A; BŠ is the greatest constant such that hAx;ÿBxi hAx; Axi for all x 2 V : 8. If A is invertible, and V ˆ R , then ÿ1 ÿ1 m ‰A; BŠˆ m ‰I ; BA Šˆ m‰ÿBA Š : V V n 9. kBxk max…m ‰A;ÿBŠ; m ‰A; BŠ†, for all x 2 V . V V 10. For regular matrices T and T, ÿ1 ÿ1 ~ ~ m …A; B†ˆ m …TA; TB†ˆ m …TAT ; TBT † V;T V TV and if V is T -invariant, ÿ1 ÿ1 ~ ~ m …A; B†ˆ m …TAT ; TBT † : V;T V Again for the euclidean norm there is a nice characterization. Proposition 3.5 Consider the n s matrix S whose columns are a basis of V . Then t t …B A‡ A B† t t t m ‰A; BŠˆ  S A AS;ÿS S : V;2 min 4 APPLICATION TO MECHANICAL SYSTEMS We consider the system 00 t u ‡ Au‡ B z ˆ f …27† Bu ˆ f …28† 2 252 I. HIGUERAS n m where u 2 R , z 2 R , A is an n n matrix such that m ‰AŠˆ > 0, i.e., Ker B hAu; ui kuk for all u 2 Ker B and B is an m n matrix such that l‰BŠˆ > 0, i.e., t m kB zk kzk for all u 2 R : As it has been pointed out in Remark 3.2, the condition imposed on B implies t t that l‰BB Š > 0 and hence the matrix BB is regular. Observe that we have dropped the condition on the continuity of A. 4.1 Perturbation Index We begin analyzing the system (28)±(28). Our aim is to obtain the perturbation index of the problem as well as to derive some information on the qualitative behavior of the solution. The DAE (27)±(28) is a particular case of the ones considered in [7]. The matrix pencil of the order 1 problem is 00 1 0 11 I 0 ÿI 0 ~ ~ @@ A @ AA A; B ˆ I ; A 0 B 0 B 00 The admissible subspace which contains the solution of the homogeneous problem is 2n‡m t ÿ1 V ˆf…u; v; z†2 R : u; v 2 Ker B ; zˆÿ…BB † BAug ; and the l.u.b. logarithmic norm was 0 ÿI ~ ~ M A; B ˆ M  I; V V A 0 0 I ˆ M ÿ A 0 ÿ1 2n t t with V ˆf…u; v†2 R : u; v 2 Ker Bg and :ˆ I ÿ B …BB † B. Observe that  is a projector onto Ker B, i.e.,  ˆ  and im  ˆ Ker B. Observe that with the euclidean inner product norm, the projector  is an orthogonal projector. STIFF INDEX-3 DAEs 253 ~ ~ A suf®cient condition for asymptotic stability of …u; v† is M ‰A; BŠ < 0, or in terms of the g.l.b. logarithmic norm, 0 ÿI m > 0 : A 0 If we have an inner product norm (Proposition 3.4-7)), 0 ÿI hAu; viÿhu; vi m ˆ min 2 2 u;v2Ker B A 0 kuk ‡kvk hAu; uiÿhu; ui min u2Ker B 2kuk 1 1 ˆ…† m ‰AŠÿ 1ˆ…† m ‰AŠÿ 1 ; Ker B Ker B 2 2 where we have used the orthogonality of . On the following we impose ~ ~ m A; B > 0 : …29† Observe that this condition implies that m ‰AŠ > 1 ; Ker B which is stronger than the Ker B-ellipticity of the map A (3). Recall that on the other hand we do not impose the continuity of A. Remark 4.1 In this part we deal with the ®nite dimensional index-3 DAE but having in mind PDAEs, i.e., _ the general operator index-3 DAE (7)±(8). Thus we must have care so that all the steps we do in the decoupling be also valid for that system. & First we are going to ®nd how the solution is decoupled in (27)±(28). It is well known that the differential index as well as the perturbation index of this DAE is 3. Recall that although in the PDAE framework the differential index has no sense, there is not any problem in de®ning the perturbation index. Thus we are going to decouple the DAE to prove that the perturbation index is 3 but using properties that are also valid for the operator DAE. This decoupling will be extremely useful for the analysis of the numerical solution. We consider consistent initial conditions, 0 0 u…t †ˆ u u…t †ˆ u 0 0 0 0 254 I. HIGUERAS satisfying 0 0 Bu ˆ f …t † ; Bu ˆ f …t † : 0 2 0 0 0 2 We begin with standard computations to obtain ÿ1 ÿ1 ÿ1 t t t 00 zˆÿ…BB † BAu‡…BB † Bf ÿ…BB † f …30† and substituting it into (27) we get ÿ1 00 t t 00 u ‡  Au ˆ f ‡ B …BB † f : …31† Any vector u can be written as u ˆ u‡…I ÿ †u. If we multiply (31) by I ÿ  and , and denote u :ˆ…I ÿ †u, u :ˆ u, we obtain 1 2 00 t t ÿ1 00 u ˆ B …BB † f 1 2 u ‡ Au ˆ f The initial condition is 0 0 Bu ˆ f …t † ; Bu ˆ f …t † ; 1;0 2 0 0 1;0 2 and thus the ®rst equation implies ÿ1 t t u ˆ B …BB † f : 1 2 From the second one we obtain ÿ1 00 t t u ‡ Au ˆ  f ÿ AB …BB † f : …32† 2 1 2 Recall that taking into account the origin of this ODE, when we deal with it, we must consider initial conditions satisfying …I ÿ †u …t †ˆ 0. 2 0 In order to study the qualitative the behavior of the solutions, it is suf®cient consider ODEs of the form (32). Thus we consider order 1 problem u ˆ v …33† v ˆÿAu‡  p …34† This is an ordinary differential equation with Im  as invariant subspace, i.e., if u…t †, v…t †2 Im , then u…t†, v…t†2 Im  for all t. To see this, we multiply 0 0 the second equation by I ÿ , to obtain …I ÿ †v ˆ 0, and hence …I ÿ † v…t†ˆ…I ÿ †v…t †. Hence if …I ÿ †v…t †ˆ 0, then …I ÿ †v…t†ˆ 0 and 0 0 therefore also …I ÿ †v …t†ˆ 0; equivalently, we obtain v…t†ˆ v…t† and STIFF INDEX-3 DAEs 255 u…t†ˆ u…t†. An important remark at this point is that the stability of (33)±(34) must be studied only on Im . This situation is not new, for example it also appears when the nonlinear index-1 DAE is decoupled [5]. Now the standard perturbations bounds for ODEs applied to (33)±(34) with invariant subspaces gives a…t† a…t† ÿa…s† k…u…t†; v…t††k  e k…u…t †; v…t †k‡ e e kp…s†k ds t  t 0 0 0 …35† ÿ1 t t where a…t†ˆ…tÿ t †m ‰A; BŠ, and p…t†ˆ f …t†ÿ  AB …BB † f …t†. 0 V 1 2 Therefore, for the solution of (27)±(28) we have kukku …t†k‡ku …t†k 2 1 a…t† 0 e…† ku…t †k‡ku…t †k 0 0 ÿ1 a…t† ÿa…s† t t ‡ e e kp…s†k ds‡kB …BB † f k ÿ1 ÿ1 ÿ1 t t t 00 kzkk…BB † BAkkuk‡k…BB † Bf ÿ…BB † f k In particular we observe that in the bounds we have f , and thus the perturbation index is 3. Condition (29) gives information on the asymp- totic stability of the problem. Recall that these bounds are also valid for (7)±(8). 4.2 Numerical Solutions by Runge±Kutta Methods In this section we show how, given (27)±(28), we can give convergence estimates independently of the stiffness of the matrix A, we will use the l.u.b. or g.l.b. constant instead. Recall that the problem we are considering is a simple lineal problem but we use techniques that can be used for for (7)±(8). We consider an s-stages Runge±Kutta method with coef®cients…b ;A†.We assume that the matrix A is regular and the method is stif¯y accurate, i.e., ÿ1 a ˆ b , i ˆ 1; .. . ; s. We denote t ˆ t ‡ c h, and ˆ…A † . As usual, si i ni nÿ1 i ij ij we assume that A l ˆ c and thus c ˆ 1. Given an approximation …u ; v ; z † of the solution at t , the new approximation …u ; v ; z † nÿ1 nÿ1 nÿ1 nÿ1 n n n at t ˆ t ‡ h is given by n nÿ1 u ˆ U ; v ˆ V ; z ˆ Z ; n ns n ns n ns 256 I. HIGUERAS where the internal stages U , V , Z , i ˆ 1; .. . ; s are obtained solving the ni ni ni system ‰UŠ ˆ V ni ni 0 t ‰VŠ ‡ AU ‡ B Z ˆ f ni ni 1i ni BU ˆ f i ˆ 1; .. . ; s ni 2i with the internal derivatives de®ned by s s X X ÿ ÿ 1 1 0 0 ‰UŠ ˆ U ÿ u ; ‰VŠ ˆ V ÿ v ij nj nÿ1 ij nj nÿ1 ni ni h h jˆ1 jˆ1 i ˆ 1; .. . ; s : We can state the following result. The proof is given in the Appendix. Theorem 4.1 Let an stif¯y accurate Runge±Kutta method …A; b † with positive weights be algebraically stable. Suppose that a positive diagonal matrix D exists such that …A† > 0. If q is the order stage of the method, then the method is convergent with order at least q for the variable u, q for the variable v, i.e., u , and qÿ 1 for the variable z . Thus for example, for the RadauIIA, we have that the orders are: s for the variable u, s for v, and sÿ 1 for z. The nonstiff order is 2sÿ 1, s and sÿ 1 respectively. Remark 4.2 Recall that Theorem 2.2, and therefore Theorem 4.1, is valid on the class F . Problems in this class must have smooth derivatives ([2, (7.2.2)]). Thus in particular it is not valid for highly oscillatory problems 00 2 like u ‡ w u ˆ 0. Roughly speaking, smoothness implies that the effective order stage is the one given by the algebraic conditions C…q†. & 5 CONCLUSIONS In this paper we have used logarithmic norms to obtain the perturbation index and convergence results for stif¯y accurate Runge±Kutta methods for a linear index-3 DAE coming from mechanical systems. This is the natural way if we want to make use of the experience of stiff ODEs [15, 2]. The decoupling for the linear index-3 DAE coming from mechanical systems is also valid for the general operator DAE (7)±(8), and thus we can de®ne the perturbation index STIFF INDEX-3 DAEs 257 also for these problems. Concerning to the numerical solution, the decoupling done is also valid for (7)±(8) but in order to ensure convergence, a detailed revision of the proofs of numerical methods for stiff ODEs [2] must be done. The open question is if properties on Kronecker products are also valid when a ®nite dimensional space and a Hilbert space are involved; these properties are implicitly used in the proofs. In any case, there is not any problem to deal with the ®nite dimensional case and, as it was expected, there is an order reduction. Taking into account that logarithmic norms are also de®ned for nonlinear maps and for unbounded operators, we conclude that this technique can also be used not only for nonlinear DAEs but also for DAEs expressed in terms of operators. This means that we have a tool to analyze solutions and numerical discretization in PDAEs, that together with the experience of theory of stiff ODEs [15, 2] may give good results. REFERENCES 1. Desoer, C. and Haneda, H.: The Measure of a Matrix as a Tool to Analyze Computer Algorithms for Circuit Analysis. IEEE Trans. Circuit Theory 19 (1972), pp. 480±486. 2. Dekker, K. and Verwer, J.G.: Stability of Runge±Kutta Methods for Stiff Nonlinear Differential Equations. North-Holland, Amsterdam, 1984. 3. Garcõ Âa-Celayeta, B.: Stability for Differential Algebraic Equations. PhD. Universidad Pu  blica de Navarra (1998). 4. Garcõ Âa-Celayeta, B. and Higueras, I.: Runge±Kutta Methods for DAEs. A New Approach. J. Comput. Appl. Math. 111 (1999), pp. 49±61. 5. Griepentrog, E. and Ma Èrz, R.: Differential Algebraic Equations and Their Numerical Treatment. Teubner Texte zur Mathematik 88, Leipzig, 1986. 6. Hairer, L., Lubich, Ch. and Roche, M.: The Numerical Solution of Differential Algebraic Equations by Runge±Kutta Methods. Lectures Notes in Mathematics 1409, Springer, Berlin, 1989. 7. Higueras, I. and Garcõa-Celayeta, B.: Logarithmic Norms for Matrix Pencils. SIAM J. Matrix Anal. Appl., 20 (1999), pp. 646±666. 8. Higueras, I. and So È derlind, G.: Logarithmic Norms and Nonlinear DAE Stability. Preprint Opto Matema Âtica e Informa Âtica. Universidad Pu  blica de Navarra Nr. 5±2001, 2001. 9. Jay, L.O.: Me Âthodes du type Runge±Kutta pour des e Âquations diffe Ârentielles alge Âbriques  Á d'index 3 avec des applications aux systemes Hamiltoniens. PhD. University of Geneve (1994). 10. Lubich, Ch.: Integration of Stiff Mechanical Systems by Runge±Kutta Methods. ZAMP,44 (1993), pp. 1022±1053. 11. Simeon, B.: Numerische Simulation gekoppelter Systeme von partiellen und differential- algebraischen Gleichungen in der Mehrk operdynamic. Habilitationsschrift, Universitat Karlsruhe, 1999. 258 I. HIGUERAS 12. So È derlind, G.: On Nonlinear Difference and Differential Equations. BIT 24 (1984), pp. 667±680. 13. So È derlind, G.: Bounds on Nonlinear Operators in Finite-dimensional Banach Spaces. Numer. Math. 50 (1986), pp. 27±44. 14. Soderlind, G.: Seminar Course Stability Concepts in Numerical Analysis. University of Lund, 1999. 15. Strom, T.: On Logarithmic Norms. SIAM J. Numer. Anal. 12 (1975), pp. 741±753. APPENDIX Proof of Theorem 4.1. In vectorial form the method can be written as ‰UŠ ˆ V …36† ‰VŠ ‡D U ‡D t Z ˆ f …37† A n B n 1 D U ˆ f i ˆ 1; .. . ; s …38† B n 2 where we have denoted U ˆ…U ; .. . ; U †, and in a similar way V , Z , n n1 ns n n 0 0 0 0 0 ‰UŠ ˆ…‰UŠ ; .. . ;‰UŠ † and in a similar way‰VŠ ,‰ZŠ ,D ˆ I A, with A s n n1 ns n n the Kronecker product [2]. In vectorial form, the internal derivatives can be expressed as ÿ1 ‰UŠ ˆ …A I†…† U ÿ l u ; n nÿ1 0 ÿ1 ‰VŠ ˆ …A I†…† V ÿ l v ; n nÿ1 with lˆ…1; .. . ; 1†2 R . In terms of the internal stages, (36)±(38) is ÿ1 …A I†…† U ÿ l u ˆ hV n nÿ1 n ÿ1 …A I†…† V ÿ l v ‡ hD U ‡ hD t Z ˆ hf n nÿ1 A n B n 1 D U ˆ f : B n 2 As usual, in order to study convergence, we consider the perturbed system ÿ1 …A I†…† U ÿ l u ˆ hV ‡ R n nÿ1 n 1 ÿ1 …A I†…† V ÿ l v ‡ hD U ‡ hD t Z ˆ hf ‡ R n nÿ1 A n B n 1 2 D U ˆ f ‡ R B n 2 3 STIFF INDEX-3 DAEs 259 We denote E ˆ U ÿ U , and in a similar way E , E , and e ˆ u ÿu , u;n n n v;n z;n u;n n n and in a similar way e , e , to obtain v;n z;n ÿ1 …A I† E ÿ l e ˆ hE ; n‡ R …39† u;n u;nÿ1 v 1 ÿ1 …A I† E ÿ l e ‡ hD E ‡ hD E ˆ R …40† v;n v;nÿ1 A u;n B z;n 2 D E ˆ R …41† B u;n 3 We multiply (40) by D , to obtain ÿ1 …A I†D E ÿD …l e † ‡ hD E ‡ hD E ˆD R B v;n B v;nÿ1 BA u;n BB z;n B 2 As BB is regular, ÿ1 hE ˆÿ…A I†D ÿ1 E ÿD ÿ1 …l e † z;n t v;n t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R t t u;n 2 …BB † BA …BB † B Multiplying (39) by D , we obtain ÿ1 …A I†D E ÿD …l e † ˆ hD E ‡D R B u;n B u;nÿ1 B v;n B 1 and hence 1 1 ÿ1 D E ˆ …A I†D E ÿD …l e † ÿ D R B v;n B u;n B u;nÿ1 B 1 h h 1 1 ÿ1 ˆ …A I† R ÿD …l e † ÿ D R 3 B u;nÿ1 B 1 h h where we have used (41). Thus, ÿ2 hE ˆÿ …A I†D ÿ1 R z;n t 3 …BB † ÿ2 ‡ …A I†D ÿ1 …l e † u;nÿ1 …BB † B ÿ1 ÿ1 ÿ1 ÿ1 ‡ …A I†D R ‡…A I†D …l e † t 1 t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R …42† t u;n t 2 …BB † BA …BB † B Finally, substituting it into (40), we get ÿ1 ÿ1 …A I†E ÿ…A I†D …l e †‡ hD E v;n  v;nÿ1 A u;n ÿ2 ÿ1 ˆD R ‡ …A I†D R 2 t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h 260 I. HIGUERAS Together with (39), we have a system of the form ÿ1 …A I†E ÿ hE ˆ T …43† u;n v;n 1;n ÿ1 hD E ‡…A I†E ˆ T …44† A u;n v;n 2;n where ÿ1 T ˆ…A I†…l e †‡ R 1;n u;nÿ1 1 ÿ1 ÿ2 T ˆ…A I†D …l e †‡D R ‡ …A I†D ÿ1 R 2;n  v;nÿ1  2 t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h Multiplying (43)±(44) by  and I ÿ , we obtain ÿ1 …A I†D E ÿ hD E ˆD T …45† u;n  v;n  1;n ÿ1 hD E ‡…A I†D E ˆD T …46† A u;n  v;n  2;n and ÿ1 …A I†D E ÿ hD E ˆD T …Iÿ† u;n …Iÿ† v;n …Iÿ† 1;n ÿ1 …A I†D E ˆD T …Iÿ† v;n …Iÿ† 2;n respectively, with ÿ1 D T ˆ…A I†D …l e †‡D R 1;n  u;nÿ1  1 ÿ1 D T ˆ…A I†D …l e †‡D R 2;n  v;nÿ1  2 and ÿ1 D T ˆ…A I†D …l e †‡D R …Iÿ† 1;n …Iÿ† u;nÿ1 …Iÿ† 1 ÿ2 D T ˆ …A I†D ÿ1 R …Iÿ† 2;n t t 3 B …BB † 1 1 ÿ2 ÿ1 ÿ …A I†D …l e †ÿ …A I†D R …Iÿ† u;nÿ1 …Iÿ† 1 h h Thus, for the I ÿ  part of the error, we have D E ˆ h…A I†D T ‡…A I†D T u;n 2;n 1;n …Iÿ† …Iÿ† …Iÿ† D E ˆ…A I†D T …Iÿ† v;n …Iÿ† 2;n STIFF INDEX-3 DAEs 261 and substituting, D T and D T we obtain …Iÿ† 1;n …Iÿ† 2;n D E ˆD ÿ1 R …47† t t …Iÿ† u;n 3 B …BB † ÿ1 D E ˆ …A I†D ÿ1 R …48† v;n t t 3 …Iÿ† B …BB † 1 1 ÿ1 ÿ …A I†D …l e †ÿ D R …49† …Iÿ† u;nÿ1 …Iÿ† 1 h h In particular, if R ˆ 0, thenD E ˆ 0, and thus…I ÿ †e ˆ 0. In this 3 u;n u; n …Iÿ† case we obtain D E ˆ 0 …50† …Iÿ† u;n D E ˆÿ D R …51† …Iÿ† v;n …Iÿ† 1 In any case, the errors in these components are not propagated in the next steps. …2† Denoting E :ˆD E for the  part of the error, we have a system of the ;n ;n form ÿ1 …2† …2† …A I†E ÿ hE ˆ  …52† 1;n u;n v;n ÿ1 …2† …2† hD E ‡…A I†E ˆ  …53† A 2;n u;n v;n with ˆD T ;  ˆD T ÿ hD E : 1;n  1;n 2;n  2;n A…Iÿ† u;n Observe that this is the numerical solution with the Runge±Kutta method of the ODE, u 0 ÿI u 0 2 2 ‡ ˆ v A 0 v 0 2 2 with Im  ˆ Ker B as invariant subspace and (29) 0 I 0 ÿI ÿM ˆ m > 0 : V V ÿ A 0  A 0 ÿ1 k ‡1 t t k ‡1 1 2 If R ˆ #…h †, and R ÿ A…I ÿ †B …BB † R ˆ #…h †, the order 1 2 3 stage we have to consider in (52)±(53) is k ˆ minfk ; k g. Thus we can apply 1 2 Theorem 2.2 of the stiff theory of ODEs to obtain B-convergence of order at least the order stage, i.e., k. 262 I. HIGUERAS q‡1 q‡1 If the RK method has order stage q, then R ˆ #…h †, R ˆ #…h †. 1 2 For the exact numerical solution, R ˆ 0. Thus for the  part of the solution we have order at least q. As for the I ÿ  part of the solution we had (50)±(51), D E ˆ 0 …Iÿ† u;n D E ˆÿ D R ; …Iÿ† v;n …Iÿ† 1 and for the error in the component z, we had (42), ÿ1 ÿ1 hE ˆ …A I†D ÿ1 R ‡…A I†D ÿ1 …l e † z;n t 1 t v;nÿ1 …BB † B …BB † B ÿ hD ÿ1 E ‡D ÿ1 R t u;n t 2 …BB † BA …BB † B

Journal

Mathematical and Computer Modelling of Dynamical Systems – Taylor & Francis

Published: Jun 1, 2001

There are no references for this article.