Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 1871590, 10 pages
https://doi.org/10.1155/2017/1871590

Research Article

Numerical Methods for a Class of Differential Algebraic Equations

Lei Ren and Yuan-Ming Wang

Department of Mathematics, Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, East China Normal University, Shanghai 200241, China

Correspondence should be addressed to Lei Ren; renleipost@163.com

Received 28 November 2016; Accepted 11 April 2017; Published 8 June 2017

Academic Editor: Fazal M. Mahomed

Copyright © 2017 Lei Ren and Yuan-Ming Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper is devoted to the study of some efficient numerical methods for the differential algebraic equations (DAEs). At first, we propose a finite algorithm to compute the Drazin inverse of the time varying DAEs. Numerical experiments are presented by Drazin inverse and Radau IIA method, which illustrate that the precision of the Drazin inverse method is higher than the Radau IIA method. Then, Drazin inverse, Radau IIA, and Padé approximation are applied to the constant coefficient DAEs, respectively. Numerical results demonstrate that the Padé approximation is powerful for solving constant coefficient DAEs.

1. Introduction

The general differential algebraic equations (DAEs)

arise in many applications, circuit analysis, control theory, chemical process simulations, singular perturbation, constrained mechanical systems of rigid bodies, and so forth [1–3]. In [4–6], the numerical solutions of constant coefficient DAEs were studied by Golub et al. The linear time varying of DAE

defined on the interval is important in understanding general DAEs. It exhibits most of the behavior found in the nonlinear case that is not already present in the constant coefficient case, yet the linearity facilitates the analysis. In spite of this, many aspects of the theory and numerical solution of (2) have only been resolved within the last few years. Some questions remain open.

In [7], one canonical form for higher index linear time varying singular systems has been presented. In [8–10], with the DAEs being tractable with indexes 2 and 3 some results were proposed. It was shown that either 2 or index 3 tractability causes the homogeneous equations with respect to (2) to provide a finite dimensional space of solutions. In [11], Song had treated general higher index time varying linear DAEs in terms of matrix pencils and investigated their basic properties and their solvability.

Above these results, the numerical solutions of DAEs seldom are involved. In this paper, we shall consider the homogeneous equations with respect to DAEs:

In this paper, we shall consider the constant coefficient and time varying DAEs with different numerical methods. In Section 2, we present some definitions which will be used in the proofs of our main theorems. In Section 3, a finite algorithm for the computation of the Drazin inverse of the time varying singular matrix is presented. Drazin inverse is applied to (3) corresponding to the difference equations; we will establish sufficient and necessary conditions for the solution of (3). In Section 4, the numerical treatment with the corresponding difference format is presented; then the numerical solutions of constant coefficient DAEs are presented by Padé approximation and the implicit Runge-Kutta method. In Section 5, some numerical examples and error estimates are proposed.

2. Preliminaries

We first introduce some definitions of the Drazin inverse of a matrix.

Definition 1. Let . If the smallest positive integer such that

holds it is called the index of and is denoted by . If A is nonsingular, then ; else if is singular, then .

Definition 2. Let and . Then the matrix satisfying

 ()

 (2)

 (3)

is called the Drazin inverse of and is denoted by .

Consider the following singular differential and corresponding difference equations

Suppose ; is sufficient small. Then

Therefore, the approximation of is given by

Let ; we can get the difference equations like (8). When is a singular matrix, things can happen that are impossible when exists.

Definition 3. For and , the vector is called a consistent initial vector associated with for the equation when the initial problem , , , processes at least one solution.

Definition 4. The equation is said to tractable at the point if the initial value problem , has a unique solution for each consistent initial vector associated with .

Definition 5. For and , the vector is called a consistent initial vector for the difference equation if the initial value problem , , , has a solution for .

Definition 6. The difference equation is said to be tractable if the initial value problem , , , has a unique solution for each consistent initial vector .

3. Time Varying Differential Algebraic Equations

Lemma 7 (see [12]). For , , and

is characteristic polynomial of , where . Then , , and

where .

Next, the algorithm for the computation of the Drazin inverse of a polynomial matrix is presented as follows.

Step 1. and are determined by the recursive relationship as follows:

and initial conditions are

Step 2. From , we can see that , . Then

The algorithm could be performed by the symbol computation package of Matlab.

We are now ready to give solutions on corresponding singular difference equation. First of all, in order to establish a sufficient and necessary condition, the following lemma is important.

Lemma 8. Let . Assume that there exists , such that exists, and let

Then .

Proof. If there exists such that exists, then

Thus
Consequently, .

Theorem 9. For , the homogeneous difference equation

is tractable if and only if there exists such that exists.

Proof. We first prove the sufficiency. Let and be defined as those in (27). Clearly is tractable if and only if is tractable.
Since , invertible matrix , so that

where is invertible and is nilpotent of index . Let . Then the differential equation becomes

or

Since is invertible, is tractable. Thus it suffices to show that the second equation of is tractable. Let and multiply (32) by . Then , and hence . So . Multiply (32) again by . Then , so . Continuing in this manner, we get and is trivially tractable.
For the necessity, assume that is tractable. We need to show that there is a such that is invertible. Assume that this is not true. Then is singular for all . This means that, for each , there is a vector such that and . Let be a finite linearly dependent set of such vectors. Let and let be such that , where not all are 0. Then is not identically zero and is easily seen to be a solution of (18). However, . Thus there are two different solutions, namely, and 0, which satisfy the initial condition . Therefore, it is not tractable at , which contradicts our hypothesis. Hence, exists for some .

The next lemma will be used to show that the solution of the difference equation is independent of the scalar which is used in the expressions and .

Lemma 10. Assume that are such that there exists such that exists. Let , . For all for which and exist, the following statements are true:

Theorem 11. If the homogeneous equation is tractable, then the general solution is given by

where and and is such that exists. Furthermore, is a consistent initial vector for if and only if , where . In this case this unique solution subject to is given by , .

Proof. Since (3) is tractable, multiplying by gives the equivalent equation . After a similarity we get, as in the proof of Theorem 9, that

The difference equation is equivalent to the pair of equations

Since is invertible, the unique solution of the second equation of (25) was . But the first equation of (25) is consistent for any and the unique solution is In terms of the original variables we have

where is arbitrary.

4. Linear Coefficient Differential Algebraic Equations

In order to establish a sufficient and necessary condition for linear coefficient DAEs, the following lemma is important.

Lemma 12 (see [13]). Let . Assume that there exists , such that exists, and let

Then .

Theorem 13. For , the homogeneous difference equation

is tractable if and only if there exists such that exists.

Proof. The proof follows from Theorem 9.

The next lemma will be used to show that the solution of the difference equation is independent of the scalar which is used in the expressions and .

Lemma 14. Assume that are such that there exists such that exists. Let , . For all for which and exist, the following statements are true:

At last, we give the general solution of the singular difference equation.

Theorem 15 (see [14]). If the homogeneous equation is tractable, then the general solution is given by

where and and is such that exists. Furthermore, is a consistent initial vector if and only if , where . In this case this unique solution subject to is given by .

Proof. The proof follows from Theorem 11

4.1. Radau IIA Methods

We will introduce the IRK which would be applied concisely. IRK methods play an important role for the numerical solution of DAEs. Due to their one-step nature, IRK methods are potentially more efficient for these problems than multistep methods because multistep methods must be restarted, usually at low order, after every discontinuity, whereas IRK methods can restart at a higher order.

The -stage implicit Runge-Kutta method applied to the general nonlinear DAE of the form (2) is defined by

The method is often denoted by the shorthand notation or Butcher diagram:

4.2. Padé Approximation [15]

A differential algebraic equation has the form

with initial values

where and are vector functions for which we assumed sufficient differentiability.

We assume the solution has the form

where is a vector function which is the same size as and . Substituting (38) into (36) and neglecting higher order term, we have the linear equation of in the form

where and are constant matrixes. Solving (39), the coefficients of can be determined. Repeating the above procedure for higher order terms, we can get the arbitrary order power series of the solutions for (2).

The power series given by above procedure can be transformed into Padé series and we have numerical solution of differential algebraic equation into (2); the specific steps can be followed by (3.1)–(3.7) in [15].

5. Numerical Experiments

Example 1. We consider the following differential algebraic equation:

where

and we take ; the analytical solution of (40) is .

(1) Drazin Inverse. Suppose ; is sufficiently small. Then

so on the point the approximate value of is given by

Now we use the difference rule to solve the equation. Set

and then

, so such that exists; we can see the homogeneous equation is tractable. For , we have such that exists. Thus

The eigenvalues of are , , and , so that could be computed by Theorem  7.5.2 in [13]:

According to Theorem 11, the general solution is given by

where , which is a consistent initial vector.

(2) Radau IIA Methods. We consider Radau IIA as follows:

The 2-stage implicit Runge-Kutta method is applied to (40):

We can rewrite (50) as follows:

where the coefficient matrix is nonsingular. So, we can compute the numerical solutions by the following equation:

(3) Padé Approximation. From initial values , the solutions of (40) can be supposed as

Substituting (53) into (40) and neglecting higher order terms, we have

The linear equation can be given in the following:

Solving this linear equation we have ; then the solutions of (40) can be supposed as

In the same manner, substituting (56) into (40) and neglecting higher order terms, then we have

Similar to (53), solving this linear equation, we have ; then the solutions of (40) can be supposed as

Repeating the above procedure we have

The power series can be transformed into the following Padé series:

If the power series converges very fast then can be eliminated in Padé series.

In Tables 1–5, exact solutions, numerical solutions, and errors are illustrated.

Table 1: Exact solution and numerical solution of .

( = 0.1)

of Drazin inv
of Radau
of Padé approx

0.1
0.3448231843
0.3443666667
0.3448231345
0.3448231842

0.2
0.3571051353
0.3563563033
0.3571050289
0.3571051354

0.3
0.3702337930
0.3691464551
0.3702336224
0.3702337931

0.4
0.3842675287
0.3827904314
0.3842672855
0.3842675289

0.5
0.3992687375
0.3973451205
0.3992684126
0.3992687378

0.6
0.4153041163
0.4128712011
0.4153036995
0.4153041166

0.7
0.4324449597
0.4294333945
0.4324444400
0.4324449593

0.8
0.4507674776
0.4471007347
0.4507668426
0.4507674772

0.9
0.4703531334
0.4659468548
0.4703523699
0.4703531335

1.0
0.4912890068
0.4860502943
0.4912881000
0.4912890064

Table 2: Exact solution and numerical solution of .

( = 0.1)

of Drazin inv
of Radau
of Padé approx

0.1
1.0919188077
1.0887666667
1.0919184095
1.0919188072

0.2
1.1901744157
1.1832248233
1.1901735645
1.1901744153

0.3
1.2952036775
1.2839732084
1.2952023127
1.2952036776

0.4
1.4074735628
1.3914309536
1.4074716175
1.4074735624

0.5
1.5274832334
1.5060450899
1.5274806342
1.5274832337

0.6
1.6557662635
1.6282924326
1.6557629295
1.6557662632

0.7
1.7928930110
1.7586815666
1.7928888531
1.7928930116

0.8
1.9394731538
1.8977549625
1.9394680743
1.9394731533

0.9
2.0961584005
2.0460912345
2.0961522922
2.0961584002

1.0
2.2636453881
2.2043075487
2.2636381331
2.2636453884

Table 3: Exact solution and numerical solution of .

( = 0.1)

of Drazin inv
of Radau
of Padé approx

0.1
−2.1493680625
−2.1443333333
−2.1493674154
−2.1493680622

0.2
−2.3090334256
−2.2981718600
−2.3090320423
−2.3090334254

0.3
−2.4797059760
−2.4622586591
−2.4797037581
−2.4797059762

0.4
−2.6621445395
−2.6372765649
−2.6621413784
−2.6621445398

0.5
−2.8571602544
−2.8239538938
−2.8571560306
−2.8571602546

0.6
−3.0656201782
−3.0230674884
−3.0656147604
−3.0656201787

0.7
−3.2884511428
−3.2354459509
−3.2884443862
−3.2884511424

0.8
−3.5266438749
−3.4619730912
−3.5266356207
−3.5266438745

0.9
−3.7812574008
−3.7035916057
−3.7812474748
−3.7812574003

1.0
−4.0534237556
−3.9613070007
−4.0534119663
−4.0534237552

Table 4: Errors .

( = 0.1)
of Drazin inverse
of Radau
of Padé approximation

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Table 5: Errors of Drazin inverse and Radau IIA.

( = 0.01)
of Drazin inverse
of Radau IIA

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Example 2 (see [16]). We consider the following differential algebraic equation:

where

and we take ; the analytical solution of (40) is .

(1) Drazin Inverse. The classical four-order Runge-Kutta method is applied to (40); we can get

where

Thus

Now we use the difference rule to solve the system of (40). Set

set ; we get , so we have such that exists. Thus

so that could be computed by finite algorithm (13)–(15):

According to Theorem 15, the general solution is given by

where .

Let ; we can compute the numerical solution of (40) by (33).

(2) Radau IIA. We consider the Radau IIA as follows:

Radau IIA is applied to (40):

We can rewrite (50) as follows:

where the coefficient matrix is nonsingular. So, we can compute the numerical solutions by the following equation:

See Table 5 for errors of the numerical results. Error is defined by .

6. Conclusion

The fundamental goal of this study has been to construct approximations to numerical solutions of linear DAEs. From Table 4, we know that the numerical solutions of Padé approximation approximate the exact solutions better than other methods, and the process of Padé approximation is easily implemented. Moreover, Drazin inverse is applied to solve the time varying differential algebraic equations. According to the obtained solutions we infer that Drazin inverse is a powerful tool for solving this kind of problems. From Table 5 we know that the precision of the Drazin inverse method is higher than the Radau IIA method, but the Drazin inverse method is implemented in more complex way than the Radau IIA method.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by Science and Technology Commission of Shanghai Municipality (STCSM) (no. 13dz2260400).

References

K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, North-Holland, New York, NY, USA, 1989.

S. L. Campbell, Singular System of Differential Equations, Pitman Publishing Ltd, San Francisco, 1980.

S. L. Campbell, Singular System of Differential Equations II, Pitman Publishing Ltd, San Francisco, 1982.

G. H. Golub and C. F. Van, Matrix Computations, Johns Hopkins University Press, Baltimore, 1983.

B. Kagstrom and A. Ruhe, Matrix Pencils, Springer-Verlag, 1983.

J. H. Wilkinson, “The practical significance of the drazin inverse,” in Recent Applications of Generalized Inverses, S. L. Campbell, Ed., pp. 82–99, Pitman, 1982.

S. L. Campbell, “One canonical form for higher-index linear time-varying singular systems,” Circuits, Systems, and Signal Processing, vol. 2, no. 3, pp. 311–326, 1983.

R. März, “Index-2 differential-algebraic equations,” Results in Mathematics, vol. 15, pp. 149–171, 1989.

R. März, “Some new results concerning index-3 differential-algebraic equations,” Journal of Mathematical Analysis and Applications, vol. 140, no. 1, pp. 177–199, 1989.

R. März, “Higher-index differential-algebraic equations: analysis and numerical treatment,” Numerical Analysis, vol. 24, pp. 199–222, 1990.

Y. Z. Sun, “Solvability of higher index time-varying linear differential-algebraic equations,” Acta Mathematica Scientia, vol. B, no. 1, pp. 77–92, 2001.

J. Gao and G. Wang, “Two algorithms of Drazin inverse of a polynomial matrix,” Joural of Shanghai Normal University, vol. 6, (Science Edition), 2002.

S. L. Campbell and C. D. Meyer, Generalized Inverse of Linear Transformations, Pitman, Great Britain, 1979.

S. L. Campbell and C. D. Meyer, Generalized Inverses of Linear Transformations, Pitman Publishing Ltd, London, 2009.

C. Ercan and B. Mustafa, “Arbitrary order numerical method for solving differential-algebraic equation by Pade series,” Applied Mathematics and Computation, vol. 137, pp. 57–65, 2003.

S. Karimi Vanani and A. Aminataei, “Numerical solution of differential algebraic equations using a multiquadric approximation scheme,” Mathematical and Computer Modelling, vol. 53, no. 5-6, pp. 659–666, 2011.