- Email: [email protected]
Index reduction for differential–algebraic equations by substitution method
Available online at www.sciencedirect.com
Linear Algebra and its Applications 429 (2008) 2268–2277 www.elsevier.com/locate/laa
Index reduction for differential–algebraic equations by substitution method Mizuyo Takamatsu a,∗ , Satoru Iwata b a Graduate School of Information Science and Technology, University of Tokyo, Tokyo 113-8656, Japan b Research Institute for Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan
Received 26 September 2007; accepted 23 June 2008 Available online 8 August 2008 Submitted by V. Mehrmann
Abstract Differential–algebraic equations (DAEs) naturally arise in many applications, but present numerical and analytical difficulties. The index of a DAE is a measure of the degree of numerical difficulty. In general, the higher the index is, the more difficult it is to solve the DAE. Therefore, it is desirable to transform the original DAE into an equivalent DAE with lower index. In this paper, we propose an index reduction method for linear DAEs with constant coefficients. The method is applicable to any DAE having at most one derivative per equality. In contrast to the other existing methods, it does not introduce any additional variables. Exploiting a combinatorial property of degrees of minors in polynomial matrices, we show that the method always reduces the index exactly by one. Thus the paper exhibits an application of combinatorial matrix theory to numerical analysis of DAEs. © 2008 Elsevier Inc. All rights reserved. AMS classification: 15A22; 34A09; 65L80 Keywords: Combinatorial matrix theory; Differential–algebraic equations; Index reduction; Matrix pencil; Kronecker canonical form
∗
Corresponding author. E-mail addresses: mizuyo− [email protected] (M. Takamatsu), [email protected] (S. Iwata).
0024-3795/$ - see front matter ( 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.laa.2008.06.025
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2269
1. Introduction Dynamical systems such as electric circuits [8,9,24], mechanical systems [22,25], and chemical plants [4,20,21,26] are often described by differential–algebraic equations (DAEs), which consist of algebraic equations and differential operations. DAEs present numerical and analytical difficulties which do not occur with ordinary differential equations (ODEs). This paper focuses on linear DAEs with constant coefficients having at most one derivative per equality. Several numerical methods have been developed for solving DAEs. For example, Gear [6] proposed the backward difference formulae (BDF), which were implemented in the DASSL code by Petzold (cf. [1]). Hairer and Wanner [10] implemented an implicit Runge–Kutta method in their RADAU5 code. The index concept plays an important role in the analysis of DAEs. The index is a measure of the degree of difficulty in the numerical solution. Roughly speaking, the higher the index is, the more difficult it is to solve the DAE. While many different concepts exist to assign an index to a DAE such as the differentiation index [1,3,10], the perturbation index [2,10], the strangeness index [14], and the tractability index [16,24], we focus on the Kronecker index in this paper. In the case of linear DAEs with constant coefficients, all these indices are equal except for the strangeness index [2,14,23]. In order to transform a DAE into an alternative form easier to solve, some index reduction methods have been developed [7,15,17]. These methods introduce additional variables, which leads to a drawback that the resulting DAE is a larger system than the original one. This paper presents a new index reduction method, called the substitution method, for linear DAEs with constant coefficients dx(t) = Ax(t) + f(t), (1) dt where E and A are constant matrices. Circuit analysis methods, such as the modified nodal analysis [11,24], are considered in the framework of this method. The substitution method is shown to reduce the index of DAEs by one without introducing any additional variables, provided that E has at most one nonzero entry in each row. This class of DAEs includes the semi-explicit form and circuit equations of most linear time-invariant circuits free from mutual inductors. The index of the linear DAE with constant coefficients (1) is given by the highest degrees of minors of specified orders in the associated polynomial matrix. The degrees of minors in polynomial matrices satisfy certain exchange property. We exploit this combinatorial property to prove our results. Thus this paper exhibits a connection between combinatorial matrix theory and analysis of DAEs. The organization of this paper is as follows. In Section 2, we explain a matrix pencil and the definition of the index. Section 3 introduces the substitution method. We describe the proposed method for index reduction in Section 4. A numerical example is given in Section 5. Section 6 concludes this paper. E
2. DAEs and matrix pencils For a polynomial p(s), we denote the degree of p(s) by deg p(s). A polynomial matrix P (s) = (pij (s)) with deg pij (s) 1 for all (i, j ) is called a matrix pencil. Obviously, a matrix pencil P (s) can be represented as P (s) = sE − A in terms of a pair of constant matrices E and A. A matrix pencil P (s) is said to be regular if P (s) is square and det P (s) is a nonvanishing polynomial.
2270
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
With the use of the Laplace transformation, the DAE in the form of (1) is expressed by the matrix pencil P (s) = sE − A as P (s)x˜ (s) = f˜ (s) + Ex(0), where s is the variable for the Laplace transform that corresponds to d/dt, the differentiation with respect to time. The linear DAE with constant coefficients (1) is solvable if and only if P (s) is a regular matrix pencil [1, Theorem 2.3.1]. The reader is referred to [1, Definition 2.2.1] for the precise definition of solvability. Thus we assume that P (s) is a regular matrix pencil throughout this paper. It is known that a regular matrix pencil can be brought into the Kronecker canonical form, which determines the index. Let Nμ denote a μ × μ matrix pencil defined by ⎛ ⎞ 1 s 0 ··· 0 ⎜ .⎟ .. ⎜0 . .. ⎟ 1 s ⎜ ⎟ ⎜ ⎟ .. .. N μ = ⎜0 . . . 0 0⎟ ⎜ ⎟ ⎜. . ⎟ .. ... ⎝ .. 1 s⎠ 0 ··· 0 0 1 Theorem 1 [5, Chapter XII, Theorem 3]. For an n × n regular matrix pencil P (s), there exist nonsingular constant matrices U and V which transform P (s) into the Kronecker canonical form: ⎛ ⎞ sIμ0 + Jμ0 O O ··· O ⎜ O Nμ1 O ··· O ⎟ ⎜ ⎟ ⎜ .. ⎟ . . ⎜ . O O Nμ2 . ⎟ U P (s)V = ⎜ ⎟, ⎜ ⎟ .. .. . . . . ⎝ . . . . O ⎠ O O ··· O Nμb where μ1 μ2 · · · μb ,
μ0 + μ1 + μ2 + · · · + μb = n,
and Jμ0 is a μ0 × μ0 Jordan canonical form. The matrices Nμi (i = 1, . . . , b) are called the nilpotent blocks. The maximum size μ1 of them is the Kronecker index, denoted by ν(P ). Obviously, ODEs have index zero, and algebraic equations have index one. For a polynomial matrix P (s) = (pij (s)), the row set and the column set are denoted by R and C, i.e., P (s) = (pij (s)|i ∈ R, j ∈ C). For I ⊆ R and J ⊆ C, P (s)[I, J ] = (pij (s)|i ∈ I, j ∈ J ) means the submatrix of P (s) with row set I and column set J . Furthermore, we denote w(I, J ) = deg det P (s)[I, J ], where w(∅, ∅) = 0 by convention. Then w enjoys the following property. Lemma 2 [19, pp. 287–289]. Let P (s) be a matrix pencil with row set R and column set C. For any (I, J ) ∈ and (I , J ) ∈ , where = {(I, J )||I | = |J |, I ⊆ R, J ⊆ C}, both (VB-1) and (VB-2) below hold: (VB-1) For any i ∈ I \ I , at least one of the following two assertions holds: (1a) There exists j ∈ J \ J such that w(I, J ) + w(I , J ) w(I \ {i}, J \ {j }) + w(I ∪ {i}, J ∪ {j }). (1b) There exists h ∈ I \ I such that
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2271
w(I, J ) + w(I , J ) w(I \ {i} ∪ {h}, J ) + w(I \ {h} ∪ {i}, J ). (VB-2) For any j ∈ J \ J , at least one of the following two assertions holds: (2a) There exists i ∈ I \ I such that w(I, J ) + w(I , J ) w(I \ {i}, J \ {j }) + w(I ∪ {i}, J ∪ {j }). (2b) There exists l ∈ J \ J such that w(I, J ) + w(I , J ) w(I, J \ {j } ∪ {l}) + w(I , J \ {l} ∪ {j }).
Let δk (P ) denote the highest degree of a minor of order k in P (s): δk (P ) = max{w(I, J )||I | = |J | = k, I ⊆ R, J ⊆ C}. I,J
We can compute δk (P ) efficiently by combinatorial relaxation algorithms [12,13,18]. The index ν(P ) is determined from δk (P ) as follows. Theorem 3 [19, Theorem 5.1.8]. Let P (s) be an n × n regular matrix pencil. The index ν(P ) is given by ν(P ) = δn−1 (P ) − δn (P ) + 1. 3. Substitution method In this section, we introduce the substitution method for solving linear DAEs with constant coefficients. The substitution method eliminates some variables by replacement to obtain a smaller system than the original one. This is familiar as the solution method for simultaneous equations. Let P (s) be an n × n regular matrix pencil with row set R and column set C, and B be a nonsingular constant submatrix of P (s) with row set X ⊂ R and column set Y ⊂ C. We transform (s) by row operations: P (s) into P
B K(s) P (s) = L(s) M(s)
I O B K(s) B K(s) (s) = →P = , (2) L(s) M(s) O D(s) −L(s)B −1 I where K(s) = P (s)[X, C \ Y ], L(s) = P (s)[R \ X, Y ], M(s) = P (s)[R \ X, C \ Y ], and D(s) = M(s) − L(s)B −1 K(s). Note that D(s) is not necessarily a matrix pencil. Consider the DAE
d Bx1 (t) + K (3) x2 (t) = f1 (t), dt
d d L x1 (t) + M x2 (t) = f2 (t). (4) dt dt By applying the transformation shown in (2), we obtain
2272
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
d Bx1 (t) = f1 (t) − K x2 (t), dt
d d D x2 (t) = f2 (t) − L B −1 f1 (t). dt dt
(5) (6)
The outline of the substitution method is as follows: Phase 1: Solve the DAE (6) for x2 (t). Phase 2: Solve the system of linear equations (5) for x1 (t). This method relies on accurate numerical computation of the transformation (2). In particular, the index of the DAE (6) depends on the resulting D dtd . In real problems, however, B is often sparse, which makes it easier to obtain the DAE (6). In the substitution method, the numerical difficulty is determined by the index ν(D) of the DAE (6). We show that ν(D) can be expressed in terms of the degrees of minors in P (s). For each i ∈ R and j ∈ C, let dij denote the degree of det P (s)[R \ {i}, C \ {j }]. Then we have (s)[R \ {i}, C \ {j }] dij = deg det P
∀i ∈ R \ X, ∀j ∈ C,
(7)
(s)[R \ {i}, C \ {j }] into P (s)[R \ {i}, C \ {j }] by row operations because we can transform P for each i ∈ R \ X and j ∈ C. The index ν(D) can be rewritten as follows. Theorem 4. For an n × n regular matrix pencil P (s), the index of D(s) is given by ν(D) = max{dij |i ∈ R \ X, j ∈ C \ Y } − δn (P ) + 1. i,j
(8)
Proof. We denote the D(s) by m. By Theorem 3, we have ν(D) = δm−1 (D) − δm (D) + 1. size of (s) = B K(s) and that B is a constant matrix. It follows from det P (s) = det P (s) Recall that P O D(s) that (s) − deg det B = deg det P (s). δm (D) = deg det D(s) = deg det P Moreover, we have δm−1 (D) = max{deg det D(s)[I, J ]||I | = |J | = m − 1} I,J
(s)[I, J ]||I | = |J | = n − 1, I ⊇ X, J ⊇ Y } − deg det B = max{deg det P I,J
= max{dij |i ∈ R \ X, j ∈ C \ Y }, i,j
where the last step is due to (7). Thus we obtain (8).
4. Index reduction Let P (s) = sE − A be an n × n regular matrix pencil such that E has at most one nonzero entry in each row. We denote the row set of P (s) by R, and the column set by C. Moreover, we assume that ν(P ) is positive. Let Y ⊆ C be the set of indices such that their column vectors in E are
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2273
zero vectors. Since P (s)[R, Y ] has full column rank by the regularity of P (s), we can find X ⊆ R such that P (s)[X, Y ] is regular. Note that because B = P (s)[X, Y ] and L = P (s)[R \ X, Y ] are (s)[R \ X, C \ Y ] is a matrix pencil. We prove that the index of D(s) constant matrices, D(s) = P is one lower than that of P (s). Lemma 5. For each i ∈ R and each j ∈ C \ Y, we have dij < δn−1 (P ). Proof. Suppose to the contrary that there exist i ∈ R and j ∈ C \ Y such that dij = δn−1 (P ). Let h be a row such that the (h, j ) entry of E is nonzero. We put (I, J ) = ({h}, {j }) and (I , J ) = (R \ {i}, C \ {j }). By (VB-2) in Lemma 2, at least one of the following two assertions holds: (2a) It follows that h = i and w({h}, {j }) + w(R \ {i}, C \ {j }) w(∅, ∅) + w(R, C). (2b) There exists l ∈ C \ {j } such that w({h}, {j }) + w(R \ {i}, C \ {j }) w({h}, {l}) + w(R \ {i}, C \ {l}). Note that w({h}, {j }) = 1 and w(R \ {i}, C \ {j }) = dij = δn−1 (P ). If (2a) holds, then it follows from w(∅, ∅) = 0 and w(R, C) = δn (P ) that 1 + δn−1 (P ) δn (P ), which implies ν(P ) 0 by Theorem 3. This contradicts ν(P ) > 0. On the other hand, if (2b) holds, we have 1 + δn−1 (P ) w({h}, {l}) + dil . Since E has at most one nonzero entry in each row, we have w({h}, {l}) = 0. Thus we obtain 1 + δn−1 (P ) dil , which contradicts the definition of δn−1 (P ). (s)[R \ X, C \ Y ] is exactly one lower than that of P (s). Theorem 6. The index of D(s) = P Proof. By Theorems 3 and 4 and Lemma 5 ν(P ) − ν(D) = δn−1 (P ) − max{dij |i ∈ R \ X, j ∈ C \ Y } > 0. i,j
We now prove ν(D) ν(P ) − 1. It follows from Lemma 5 that there exist i ∈ R and j ∈ Y such that dij = δn−1 (P ). Suppose that there exist i ∈ R \ X and j ∈ Y such that dij = δn−1 (P ). By applying (VB-2) in Lemma 2 to (X, Y ) and (R \ {i}, C \ {j }), there exists l ∈ C \ Y such that w(X, Y ) + w(R \ {i}, C \ {j }) w(X, Y \ {j } ∪ {l}) + w(R \ {i}, C \ {l}). Note that w(X, Y ) = 0, because P (s)[R, Y ] is a constant matrix. Since P (s) is a matrix pencil and P (s)[X, Y ] is a constant matrix, w(X, Y \ {j } ∪ {l}) 1. Therefore, we have dij dil + 1, which implies ν(D) dil − δn (P ) + 1 dij − δn (P ) = ν(P ) − 1 by Theorems 3 and 4. We now consider the other case, which means that there exist i ∈ X and j ∈ Y such that dij = δn−1 (P ), and duv < δn−1 (P ) for any u ∈ R \ X and v ∈ Y . By applying (VB-1) in Lemma 2 to (X, Y ) and (R \ {i}, C \ {j }), at least one of the following assertions holds: (1a) It follows that w(X, Y ) + w(R \ {i}, C \ {j }) w(X \ {i}, Y \ {j }) + w(R, C). (1b) There exists h ∈ R \ X such that w(X, Y ) + w(R \ {i}, C \ {j }) w(X \ {i} ∪ {h}, Y ) + w(R \ {h}, C \ {j }). Since P (s)[R, Y ] is a constant matrix, we have w(X, Y ) = w(X \ {i}, Y \ {j }) = w(X \ {i} ∪ {h}, Y ) = 0.
2274
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
C
V
a
I
L
Fig. 1. Linear circuit described by circuit equations with index three.
If (1a) holds, then we have dij δn (P ). Therefore, ν(P ) = dij − δn (P ) + 1 1 by Theorem 3. It follows from the nonnegativity of ν(D) that ν(D) ν(P ) − 1. On the other hand, if (1b) holds, we have dij dhj . This contradicts the assumption that duv < δn−1 (P ) for any u ∈ R \ X and v ∈ Y . Theorem 6 implies that the index of D(s) is the same for any X with P (s)[X, Y ] being a nonsingular constant matrix. 5. Numerical example In this section, we demonstrate the proposed method in a numerical example. We use RADAU5 [10] in Matlab as the DAE solver. RADAU5 is an implementation of a fifth order implicit Runge– Kutta method with three stages (RADAU IIA). This is applicable to ODEs and DAEs with index at most three. Example 7 (Electric circuit with index three [9]). Consider a circuit with a current-controlled current source I depicted in Fig. 1. The current through I is controlled by the current through V . This circuit is described by the circuit equations with index three:
Fig. 2. The current through the inductance: numerical solutions of the original DAE (dash-dotted line), the substitution method (solid line), and the exact solution (dotted line).
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2275
Fig. 3. The error in the current through the inductance: the original DAE (dash-dotted line) and the substitution method (solid line).
c1 c2 r1 1 1 r2 ⎜ 0 0 ⎜ 0 r3 ⎜ ⎜0 0 0 r4 ⎜ ⎜ 0 0 r5 ⎜ ⎜ r6 ⎜ ⎜0 −1 0 r7 ⎝ 0 0 r8 a ⎛
c3 c 4 0 0 1 1 0 0 0 0 0 0 0 0 0 sL 1 0
c5 c6 0 0 0 0 1 −1 0 0 1 0 0 sC 0 0 0 0
c7 c8 ⎞ ⎞ ⎛˜ ⎞ ⎛ iV 0 0 0 ⎜˜ ⎟ ⎜ ⎟ 0 0⎟ ⎟ ⎜ iC ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎟ ⎜ ˜ 0 0 ⎟ ⎜ iI ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎟ ⎜ ⎟ 1 −1⎟ ⎟ ⎜ i˜L ⎟ = ⎜ 0 ⎟ . ⎟ ⎜ ⎟ ⎜ ˜ 0 0 ⎟ ⎜v˜V ⎟ ⎜V (s)⎟ ⎟ ⎜ ⎟ ⎜ ⎟ 0 0⎟ ⎟ ⎜v˜C ⎟ ⎜ 0 ⎟ ⎠ ⎝ ⎝ ⎠ 0 −1 0 ⎠ v˜I 0 0 0 v˜L
The modified nodal analysis results in a DAE with index three [9]. However, our method finds X = {r1 , r4 , r5 , r6 , r7 , r8 } and we obtain
1 D(s) = 0
and
Y = {c1 , c2 , c3 , c5 , c7 , c8 },
saC , −1
which has index two. Setting C = 5 × 10−6 , L = 8 × 10−3 , a = 0.99, and V (t) = 10 sin(200t), we numerically solve both the original and the resulting DAEs. Fig. 2 presents these two numerical solutions and the exact solution, which can be obtained analytically. In Fig. 2, the exact solution coincides with the solution of the substitution method. Fig. 3 shows the discrepancy of the two numerical solutions from the exact solution. It is observed that the index reduction effectively improves the accuracy of the numerical solution. 6. Conclusion For linear DAEs with constant coefficients, we have proposed a new index reduction method. This method is applicable to all DAEs with at most one derivative per equality, and always reduces
2276
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
the index by one without introducing any additional variables. This class of DAEs includes the semi-explicit form and circuit equations of most linear time-invariant circuits.
Acknowledgements We are grateful to Taketomo Mitsui, Kazuo Murota, and Caren Tischendorf for their helpful comments and suggestions.
References [1] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential–Algebraic Equations, second ed., SIAM, Philadelphia, 1996. [2] P. Bujakiewicz, Maximum weighted matching for high index differential algebraic equations, Doctor’s dissertation, Delft University of Technology, 1994. [3] S.L. Campbell, C.W. Gear, The index of general nonlinear DAEs, Numer. Math. 72 (1995) 173–196. [4] R. Gani, I.T. Cameron, Modelling for dynamic simulation of chemical processes: the index problem, Chem. Eng. Sci. 47 (1992) 1311–1315. [5] F.R. Gantmacher, The Theory of Matrices, Chelsea, New York, 1959. [6] C.W. Gear, Simultaneous numerical solution of differential–algebraic equations, IEEE Trans. Circuit Theory 18 (1971) 89–95. [7] C.W. Gear, Differential–algebraic equation index transformations, SIAM J. Sci. Stat. Comput. 9 (1988) 39– 47. [8] M. Günther, U. Feldmann, The DAE-index in electric circuit simulation, Math. Comput. Simul. 39 (1995) 573– 582. [9] M. Günther, P. Rentrop, The differential–algebraic index concept in electric circuit simulation, Z. Angew. Math. Mech. 76 (Suppl. 1) (1996) 91–94. [10] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, second ed., Springer-Verlag, Berlin, 1996. [11] C.W. Ho, A.E. Ruehli, P.A. Brennan, The modified nodal approach to network analysis,IEEE Trans. Circuits Systems CAS-22 (1975) 504–509. [12] S. Iwata, Computing the maximum degree of minors in matrix pencils via combinatorial relaxation, Algorithmica 36 (2003) 331–341. [13] S. Iwata, K. Murota, I. Sakuta, Primal–dual combinatorial relaxation algorithms for the maximum degree of subdeterminants, SIAM J. Sci. Comput. 17 (1996) 993–1012. [14] P. Kunkel, V. Mehrmann, Canonical forms for linear differential–algebraic equations with variable coefficients, J. Comput. Appl. Math. 56 (1994) 225–251. [15] P. Kunkel, V. Mehrmann, Index reduction for differential–algebraic equations by minimal extension, Z. Angew. Math. Mech. 84 (2004) 579–597. [16] R. März, Numerical methods for differential–algebraic equations, Acta Numer. 1 (1992) 141–198. [17] S.E. Mattsson, G. Söderlind, Index reduction in differential–algebraic equations using dummy derivatives, SIAM J. Sci. Comput. 14 (1993) 677–692. [18] K. Murota, Combinatorial relaxation algorithm for the maximum degree of subdeterminants: computing Smith–McMillan form at infinity structural indices in Kronecker form, Appl. Algebra Engrg. Comm. Comput. 6 (1995) 251–273. [19] K. Murota, Matrices and Matroids for Systems Analysis, Springer-Verlag, Berlin, 2000. [20] C.C. Pantelides, D. Gritsis, K.R. Morison, R.W.H. Sargent, The mathematical modelling of transient systems using differential–algebraic equations, Comput. Chem. Eng. 12 (1988) 449–454. [21] J.W. Ponton, P.J. Gawthrop, Systematic construction of dynamic models for phase equilibrium processes, Comput. Chem. Eng. 15 (1991) 803–808. [22] P. Rabier, W.C. Rheinboldt, Nonholonomic Motion of Rigid Mechanical Systems from a DAE Viewpoint, SIAM, Philadelphia, 2000. [23] S. Schulz, Four lectures on differential–algebraic equations, Technical Report 497, The University of Auckland, New Zealand, 2003.
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2277
[24] D.E. Schwarz, C. Tischendorf, Structural analysis of electric circuits and consequences for MNA, Internat. J. Circuit Theory Appl. 28 (2000) 131–162. [25] B. Simeon, C. Fuhrer, P. Rentrop, Differential–algebraic equations in vehicle system dynamics, Surveys Math. Indust. 1 (1991) 1–37. [26] J. Unger, A. Kröner, W. Marquardt, Structural analysis of differential–algebraic equation systems – theory and applications, Comput. Chem. Eng. 19 (1995) 867–882.
Linear Algebra and its Applications 429 (2008) 2268–2277 www.elsevier.com/locate/laa
Index reduction for differential–algebraic equations by substitution method Mizuyo Takamatsu a,∗ , Satoru Iwata b a Graduate School of Information Science and Technology, University of Tokyo, Tokyo 113-8656, Japan b Research Institute for Mathematical Sciences, Kyoto University, Kyoto 606-8502, Japan
Received 26 September 2007; accepted 23 June 2008 Available online 8 August 2008 Submitted by V. Mehrmann
Abstract Differential–algebraic equations (DAEs) naturally arise in many applications, but present numerical and analytical difficulties. The index of a DAE is a measure of the degree of numerical difficulty. In general, the higher the index is, the more difficult it is to solve the DAE. Therefore, it is desirable to transform the original DAE into an equivalent DAE with lower index. In this paper, we propose an index reduction method for linear DAEs with constant coefficients. The method is applicable to any DAE having at most one derivative per equality. In contrast to the other existing methods, it does not introduce any additional variables. Exploiting a combinatorial property of degrees of minors in polynomial matrices, we show that the method always reduces the index exactly by one. Thus the paper exhibits an application of combinatorial matrix theory to numerical analysis of DAEs. © 2008 Elsevier Inc. All rights reserved. AMS classification: 15A22; 34A09; 65L80 Keywords: Combinatorial matrix theory; Differential–algebraic equations; Index reduction; Matrix pencil; Kronecker canonical form
∗
Corresponding author. E-mail addresses: mizuyo− [email protected] (M. Takamatsu), [email protected] (S. Iwata).
0024-3795/$ - see front matter ( 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.laa.2008.06.025
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2269
1. Introduction Dynamical systems such as electric circuits [8,9,24], mechanical systems [22,25], and chemical plants [4,20,21,26] are often described by differential–algebraic equations (DAEs), which consist of algebraic equations and differential operations. DAEs present numerical and analytical difficulties which do not occur with ordinary differential equations (ODEs). This paper focuses on linear DAEs with constant coefficients having at most one derivative per equality. Several numerical methods have been developed for solving DAEs. For example, Gear [6] proposed the backward difference formulae (BDF), which were implemented in the DASSL code by Petzold (cf. [1]). Hairer and Wanner [10] implemented an implicit Runge–Kutta method in their RADAU5 code. The index concept plays an important role in the analysis of DAEs. The index is a measure of the degree of difficulty in the numerical solution. Roughly speaking, the higher the index is, the more difficult it is to solve the DAE. While many different concepts exist to assign an index to a DAE such as the differentiation index [1,3,10], the perturbation index [2,10], the strangeness index [14], and the tractability index [16,24], we focus on the Kronecker index in this paper. In the case of linear DAEs with constant coefficients, all these indices are equal except for the strangeness index [2,14,23]. In order to transform a DAE into an alternative form easier to solve, some index reduction methods have been developed [7,15,17]. These methods introduce additional variables, which leads to a drawback that the resulting DAE is a larger system than the original one. This paper presents a new index reduction method, called the substitution method, for linear DAEs with constant coefficients dx(t) = Ax(t) + f(t), (1) dt where E and A are constant matrices. Circuit analysis methods, such as the modified nodal analysis [11,24], are considered in the framework of this method. The substitution method is shown to reduce the index of DAEs by one without introducing any additional variables, provided that E has at most one nonzero entry in each row. This class of DAEs includes the semi-explicit form and circuit equations of most linear time-invariant circuits free from mutual inductors. The index of the linear DAE with constant coefficients (1) is given by the highest degrees of minors of specified orders in the associated polynomial matrix. The degrees of minors in polynomial matrices satisfy certain exchange property. We exploit this combinatorial property to prove our results. Thus this paper exhibits a connection between combinatorial matrix theory and analysis of DAEs. The organization of this paper is as follows. In Section 2, we explain a matrix pencil and the definition of the index. Section 3 introduces the substitution method. We describe the proposed method for index reduction in Section 4. A numerical example is given in Section 5. Section 6 concludes this paper. E
2. DAEs and matrix pencils For a polynomial p(s), we denote the degree of p(s) by deg p(s). A polynomial matrix P (s) = (pij (s)) with deg pij (s) 1 for all (i, j ) is called a matrix pencil. Obviously, a matrix pencil P (s) can be represented as P (s) = sE − A in terms of a pair of constant matrices E and A. A matrix pencil P (s) is said to be regular if P (s) is square and det P (s) is a nonvanishing polynomial.
2270
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
With the use of the Laplace transformation, the DAE in the form of (1) is expressed by the matrix pencil P (s) = sE − A as P (s)x˜ (s) = f˜ (s) + Ex(0), where s is the variable for the Laplace transform that corresponds to d/dt, the differentiation with respect to time. The linear DAE with constant coefficients (1) is solvable if and only if P (s) is a regular matrix pencil [1, Theorem 2.3.1]. The reader is referred to [1, Definition 2.2.1] for the precise definition of solvability. Thus we assume that P (s) is a regular matrix pencil throughout this paper. It is known that a regular matrix pencil can be brought into the Kronecker canonical form, which determines the index. Let Nμ denote a μ × μ matrix pencil defined by ⎛ ⎞ 1 s 0 ··· 0 ⎜ .⎟ .. ⎜0 . .. ⎟ 1 s ⎜ ⎟ ⎜ ⎟ .. .. N μ = ⎜0 . . . 0 0⎟ ⎜ ⎟ ⎜. . ⎟ .. ... ⎝ .. 1 s⎠ 0 ··· 0 0 1 Theorem 1 [5, Chapter XII, Theorem 3]. For an n × n regular matrix pencil P (s), there exist nonsingular constant matrices U and V which transform P (s) into the Kronecker canonical form: ⎛ ⎞ sIμ0 + Jμ0 O O ··· O ⎜ O Nμ1 O ··· O ⎟ ⎜ ⎟ ⎜ .. ⎟ . . ⎜ . O O Nμ2 . ⎟ U P (s)V = ⎜ ⎟, ⎜ ⎟ .. .. . . . . ⎝ . . . . O ⎠ O O ··· O Nμb where μ1 μ2 · · · μb ,
μ0 + μ1 + μ2 + · · · + μb = n,
and Jμ0 is a μ0 × μ0 Jordan canonical form. The matrices Nμi (i = 1, . . . , b) are called the nilpotent blocks. The maximum size μ1 of them is the Kronecker index, denoted by ν(P ). Obviously, ODEs have index zero, and algebraic equations have index one. For a polynomial matrix P (s) = (pij (s)), the row set and the column set are denoted by R and C, i.e., P (s) = (pij (s)|i ∈ R, j ∈ C). For I ⊆ R and J ⊆ C, P (s)[I, J ] = (pij (s)|i ∈ I, j ∈ J ) means the submatrix of P (s) with row set I and column set J . Furthermore, we denote w(I, J ) = deg det P (s)[I, J ], where w(∅, ∅) = 0 by convention. Then w enjoys the following property. Lemma 2 [19, pp. 287–289]. Let P (s) be a matrix pencil with row set R and column set C. For any (I, J ) ∈ and (I , J ) ∈ , where = {(I, J )||I | = |J |, I ⊆ R, J ⊆ C}, both (VB-1) and (VB-2) below hold: (VB-1) For any i ∈ I \ I , at least one of the following two assertions holds: (1a) There exists j ∈ J \ J such that w(I, J ) + w(I , J ) w(I \ {i}, J \ {j }) + w(I ∪ {i}, J ∪ {j }). (1b) There exists h ∈ I \ I such that
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2271
w(I, J ) + w(I , J ) w(I \ {i} ∪ {h}, J ) + w(I \ {h} ∪ {i}, J ). (VB-2) For any j ∈ J \ J , at least one of the following two assertions holds: (2a) There exists i ∈ I \ I such that w(I, J ) + w(I , J ) w(I \ {i}, J \ {j }) + w(I ∪ {i}, J ∪ {j }). (2b) There exists l ∈ J \ J such that w(I, J ) + w(I , J ) w(I, J \ {j } ∪ {l}) + w(I , J \ {l} ∪ {j }).
Let δk (P ) denote the highest degree of a minor of order k in P (s): δk (P ) = max{w(I, J )||I | = |J | = k, I ⊆ R, J ⊆ C}. I,J
We can compute δk (P ) efficiently by combinatorial relaxation algorithms [12,13,18]. The index ν(P ) is determined from δk (P ) as follows. Theorem 3 [19, Theorem 5.1.8]. Let P (s) be an n × n regular matrix pencil. The index ν(P ) is given by ν(P ) = δn−1 (P ) − δn (P ) + 1. 3. Substitution method In this section, we introduce the substitution method for solving linear DAEs with constant coefficients. The substitution method eliminates some variables by replacement to obtain a smaller system than the original one. This is familiar as the solution method for simultaneous equations. Let P (s) be an n × n regular matrix pencil with row set R and column set C, and B be a nonsingular constant submatrix of P (s) with row set X ⊂ R and column set Y ⊂ C. We transform (s) by row operations: P (s) into P
B K(s) P (s) = L(s) M(s)
I O B K(s) B K(s) (s) = →P = , (2) L(s) M(s) O D(s) −L(s)B −1 I where K(s) = P (s)[X, C \ Y ], L(s) = P (s)[R \ X, Y ], M(s) = P (s)[R \ X, C \ Y ], and D(s) = M(s) − L(s)B −1 K(s). Note that D(s) is not necessarily a matrix pencil. Consider the DAE
d Bx1 (t) + K (3) x2 (t) = f1 (t), dt
d d L x1 (t) + M x2 (t) = f2 (t). (4) dt dt By applying the transformation shown in (2), we obtain
2272
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
d Bx1 (t) = f1 (t) − K x2 (t), dt
d d D x2 (t) = f2 (t) − L B −1 f1 (t). dt dt
(5) (6)
The outline of the substitution method is as follows: Phase 1: Solve the DAE (6) for x2 (t). Phase 2: Solve the system of linear equations (5) for x1 (t). This method relies on accurate numerical computation of the transformation (2). In particular, the index of the DAE (6) depends on the resulting D dtd . In real problems, however, B is often sparse, which makes it easier to obtain the DAE (6). In the substitution method, the numerical difficulty is determined by the index ν(D) of the DAE (6). We show that ν(D) can be expressed in terms of the degrees of minors in P (s). For each i ∈ R and j ∈ C, let dij denote the degree of det P (s)[R \ {i}, C \ {j }]. Then we have (s)[R \ {i}, C \ {j }] dij = deg det P
∀i ∈ R \ X, ∀j ∈ C,
(7)
(s)[R \ {i}, C \ {j }] into P (s)[R \ {i}, C \ {j }] by row operations because we can transform P for each i ∈ R \ X and j ∈ C. The index ν(D) can be rewritten as follows. Theorem 4. For an n × n regular matrix pencil P (s), the index of D(s) is given by ν(D) = max{dij |i ∈ R \ X, j ∈ C \ Y } − δn (P ) + 1. i,j
(8)
Proof. We denote the D(s) by m. By Theorem 3, we have ν(D) = δm−1 (D) − δm (D) + 1. size of (s) = B K(s) and that B is a constant matrix. It follows from det P (s) = det P (s) Recall that P O D(s) that (s) − deg det B = deg det P (s). δm (D) = deg det D(s) = deg det P Moreover, we have δm−1 (D) = max{deg det D(s)[I, J ]||I | = |J | = m − 1} I,J
(s)[I, J ]||I | = |J | = n − 1, I ⊇ X, J ⊇ Y } − deg det B = max{deg det P I,J
= max{dij |i ∈ R \ X, j ∈ C \ Y }, i,j
where the last step is due to (7). Thus we obtain (8).
4. Index reduction Let P (s) = sE − A be an n × n regular matrix pencil such that E has at most one nonzero entry in each row. We denote the row set of P (s) by R, and the column set by C. Moreover, we assume that ν(P ) is positive. Let Y ⊆ C be the set of indices such that their column vectors in E are
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2273
zero vectors. Since P (s)[R, Y ] has full column rank by the regularity of P (s), we can find X ⊆ R such that P (s)[X, Y ] is regular. Note that because B = P (s)[X, Y ] and L = P (s)[R \ X, Y ] are (s)[R \ X, C \ Y ] is a matrix pencil. We prove that the index of D(s) constant matrices, D(s) = P is one lower than that of P (s). Lemma 5. For each i ∈ R and each j ∈ C \ Y, we have dij < δn−1 (P ). Proof. Suppose to the contrary that there exist i ∈ R and j ∈ C \ Y such that dij = δn−1 (P ). Let h be a row such that the (h, j ) entry of E is nonzero. We put (I, J ) = ({h}, {j }) and (I , J ) = (R \ {i}, C \ {j }). By (VB-2) in Lemma 2, at least one of the following two assertions holds: (2a) It follows that h = i and w({h}, {j }) + w(R \ {i}, C \ {j }) w(∅, ∅) + w(R, C). (2b) There exists l ∈ C \ {j } such that w({h}, {j }) + w(R \ {i}, C \ {j }) w({h}, {l}) + w(R \ {i}, C \ {l}). Note that w({h}, {j }) = 1 and w(R \ {i}, C \ {j }) = dij = δn−1 (P ). If (2a) holds, then it follows from w(∅, ∅) = 0 and w(R, C) = δn (P ) that 1 + δn−1 (P ) δn (P ), which implies ν(P ) 0 by Theorem 3. This contradicts ν(P ) > 0. On the other hand, if (2b) holds, we have 1 + δn−1 (P ) w({h}, {l}) + dil . Since E has at most one nonzero entry in each row, we have w({h}, {l}) = 0. Thus we obtain 1 + δn−1 (P ) dil , which contradicts the definition of δn−1 (P ). (s)[R \ X, C \ Y ] is exactly one lower than that of P (s). Theorem 6. The index of D(s) = P Proof. By Theorems 3 and 4 and Lemma 5 ν(P ) − ν(D) = δn−1 (P ) − max{dij |i ∈ R \ X, j ∈ C \ Y } > 0. i,j
We now prove ν(D) ν(P ) − 1. It follows from Lemma 5 that there exist i ∈ R and j ∈ Y such that dij = δn−1 (P ). Suppose that there exist i ∈ R \ X and j ∈ Y such that dij = δn−1 (P ). By applying (VB-2) in Lemma 2 to (X, Y ) and (R \ {i}, C \ {j }), there exists l ∈ C \ Y such that w(X, Y ) + w(R \ {i}, C \ {j }) w(X, Y \ {j } ∪ {l}) + w(R \ {i}, C \ {l}). Note that w(X, Y ) = 0, because P (s)[R, Y ] is a constant matrix. Since P (s) is a matrix pencil and P (s)[X, Y ] is a constant matrix, w(X, Y \ {j } ∪ {l}) 1. Therefore, we have dij dil + 1, which implies ν(D) dil − δn (P ) + 1 dij − δn (P ) = ν(P ) − 1 by Theorems 3 and 4. We now consider the other case, which means that there exist i ∈ X and j ∈ Y such that dij = δn−1 (P ), and duv < δn−1 (P ) for any u ∈ R \ X and v ∈ Y . By applying (VB-1) in Lemma 2 to (X, Y ) and (R \ {i}, C \ {j }), at least one of the following assertions holds: (1a) It follows that w(X, Y ) + w(R \ {i}, C \ {j }) w(X \ {i}, Y \ {j }) + w(R, C). (1b) There exists h ∈ R \ X such that w(X, Y ) + w(R \ {i}, C \ {j }) w(X \ {i} ∪ {h}, Y ) + w(R \ {h}, C \ {j }). Since P (s)[R, Y ] is a constant matrix, we have w(X, Y ) = w(X \ {i}, Y \ {j }) = w(X \ {i} ∪ {h}, Y ) = 0.
2274
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
C
V
a
I
L
Fig. 1. Linear circuit described by circuit equations with index three.
If (1a) holds, then we have dij δn (P ). Therefore, ν(P ) = dij − δn (P ) + 1 1 by Theorem 3. It follows from the nonnegativity of ν(D) that ν(D) ν(P ) − 1. On the other hand, if (1b) holds, we have dij dhj . This contradicts the assumption that duv < δn−1 (P ) for any u ∈ R \ X and v ∈ Y . Theorem 6 implies that the index of D(s) is the same for any X with P (s)[X, Y ] being a nonsingular constant matrix. 5. Numerical example In this section, we demonstrate the proposed method in a numerical example. We use RADAU5 [10] in Matlab as the DAE solver. RADAU5 is an implementation of a fifth order implicit Runge– Kutta method with three stages (RADAU IIA). This is applicable to ODEs and DAEs with index at most three. Example 7 (Electric circuit with index three [9]). Consider a circuit with a current-controlled current source I depicted in Fig. 1. The current through I is controlled by the current through V . This circuit is described by the circuit equations with index three:
Fig. 2. The current through the inductance: numerical solutions of the original DAE (dash-dotted line), the substitution method (solid line), and the exact solution (dotted line).
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2275
Fig. 3. The error in the current through the inductance: the original DAE (dash-dotted line) and the substitution method (solid line).
c1 c2 r1 1 1 r2 ⎜ 0 0 ⎜ 0 r3 ⎜ ⎜0 0 0 r4 ⎜ ⎜ 0 0 r5 ⎜ ⎜ r6 ⎜ ⎜0 −1 0 r7 ⎝ 0 0 r8 a ⎛
c3 c 4 0 0 1 1 0 0 0 0 0 0 0 0 0 sL 1 0
c5 c6 0 0 0 0 1 −1 0 0 1 0 0 sC 0 0 0 0
c7 c8 ⎞ ⎞ ⎛˜ ⎞ ⎛ iV 0 0 0 ⎜˜ ⎟ ⎜ ⎟ 0 0⎟ ⎟ ⎜ iC ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎟ ⎜ ˜ 0 0 ⎟ ⎜ iI ⎟ ⎜ 0 ⎟ ⎟ ⎜ ⎟ ⎜ ⎟ 1 −1⎟ ⎟ ⎜ i˜L ⎟ = ⎜ 0 ⎟ . ⎟ ⎜ ⎟ ⎜ ˜ 0 0 ⎟ ⎜v˜V ⎟ ⎜V (s)⎟ ⎟ ⎜ ⎟ ⎜ ⎟ 0 0⎟ ⎟ ⎜v˜C ⎟ ⎜ 0 ⎟ ⎠ ⎝ ⎝ ⎠ 0 −1 0 ⎠ v˜I 0 0 0 v˜L
The modified nodal analysis results in a DAE with index three [9]. However, our method finds X = {r1 , r4 , r5 , r6 , r7 , r8 } and we obtain
1 D(s) = 0
and
Y = {c1 , c2 , c3 , c5 , c7 , c8 },
saC , −1
which has index two. Setting C = 5 × 10−6 , L = 8 × 10−3 , a = 0.99, and V (t) = 10 sin(200t), we numerically solve both the original and the resulting DAEs. Fig. 2 presents these two numerical solutions and the exact solution, which can be obtained analytically. In Fig. 2, the exact solution coincides with the solution of the substitution method. Fig. 3 shows the discrepancy of the two numerical solutions from the exact solution. It is observed that the index reduction effectively improves the accuracy of the numerical solution. 6. Conclusion For linear DAEs with constant coefficients, we have proposed a new index reduction method. This method is applicable to all DAEs with at most one derivative per equality, and always reduces
2276
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
the index by one without introducing any additional variables. This class of DAEs includes the semi-explicit form and circuit equations of most linear time-invariant circuits.
Acknowledgements We are grateful to Taketomo Mitsui, Kazuo Murota, and Caren Tischendorf for their helpful comments and suggestions.
References [1] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential–Algebraic Equations, second ed., SIAM, Philadelphia, 1996. [2] P. Bujakiewicz, Maximum weighted matching for high index differential algebraic equations, Doctor’s dissertation, Delft University of Technology, 1994. [3] S.L. Campbell, C.W. Gear, The index of general nonlinear DAEs, Numer. Math. 72 (1995) 173–196. [4] R. Gani, I.T. Cameron, Modelling for dynamic simulation of chemical processes: the index problem, Chem. Eng. Sci. 47 (1992) 1311–1315. [5] F.R. Gantmacher, The Theory of Matrices, Chelsea, New York, 1959. [6] C.W. Gear, Simultaneous numerical solution of differential–algebraic equations, IEEE Trans. Circuit Theory 18 (1971) 89–95. [7] C.W. Gear, Differential–algebraic equation index transformations, SIAM J. Sci. Stat. Comput. 9 (1988) 39– 47. [8] M. Günther, U. Feldmann, The DAE-index in electric circuit simulation, Math. Comput. Simul. 39 (1995) 573– 582. [9] M. Günther, P. Rentrop, The differential–algebraic index concept in electric circuit simulation, Z. Angew. Math. Mech. 76 (Suppl. 1) (1996) 91–94. [10] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, second ed., Springer-Verlag, Berlin, 1996. [11] C.W. Ho, A.E. Ruehli, P.A. Brennan, The modified nodal approach to network analysis,IEEE Trans. Circuits Systems CAS-22 (1975) 504–509. [12] S. Iwata, Computing the maximum degree of minors in matrix pencils via combinatorial relaxation, Algorithmica 36 (2003) 331–341. [13] S. Iwata, K. Murota, I. Sakuta, Primal–dual combinatorial relaxation algorithms for the maximum degree of subdeterminants, SIAM J. Sci. Comput. 17 (1996) 993–1012. [14] P. Kunkel, V. Mehrmann, Canonical forms for linear differential–algebraic equations with variable coefficients, J. Comput. Appl. Math. 56 (1994) 225–251. [15] P. Kunkel, V. Mehrmann, Index reduction for differential–algebraic equations by minimal extension, Z. Angew. Math. Mech. 84 (2004) 579–597. [16] R. März, Numerical methods for differential–algebraic equations, Acta Numer. 1 (1992) 141–198. [17] S.E. Mattsson, G. Söderlind, Index reduction in differential–algebraic equations using dummy derivatives, SIAM J. Sci. Comput. 14 (1993) 677–692. [18] K. Murota, Combinatorial relaxation algorithm for the maximum degree of subdeterminants: computing Smith–McMillan form at infinity structural indices in Kronecker form, Appl. Algebra Engrg. Comm. Comput. 6 (1995) 251–273. [19] K. Murota, Matrices and Matroids for Systems Analysis, Springer-Verlag, Berlin, 2000. [20] C.C. Pantelides, D. Gritsis, K.R. Morison, R.W.H. Sargent, The mathematical modelling of transient systems using differential–algebraic equations, Comput. Chem. Eng. 12 (1988) 449–454. [21] J.W. Ponton, P.J. Gawthrop, Systematic construction of dynamic models for phase equilibrium processes, Comput. Chem. Eng. 15 (1991) 803–808. [22] P. Rabier, W.C. Rheinboldt, Nonholonomic Motion of Rigid Mechanical Systems from a DAE Viewpoint, SIAM, Philadelphia, 2000. [23] S. Schulz, Four lectures on differential–algebraic equations, Technical Report 497, The University of Auckland, New Zealand, 2003.
M. Takamatsu, S. Iwata / Linear Algebra and its Applications 429 (2008) 2268–2277
2277
[24] D.E. Schwarz, C. Tischendorf, Structural analysis of electric circuits and consequences for MNA, Internat. J. Circuit Theory Appl. 28 (2000) 131–162. [25] B. Simeon, C. Fuhrer, P. Rentrop, Differential–algebraic equations in vehicle system dynamics, Surveys Math. Indust. 1 (1991) 1–37. [26] J. Unger, A. Kröner, W. Marquardt, Structural analysis of differential–algebraic equation systems – theory and applications, Comput. Chem. Eng. 19 (1995) 867–882.