# Numerical Methods for a Class of Differential Algebraic Equations

Authors: Lei Ren; Yuan-Ming Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 1. Introduction

The general differential algebraic equations (DAEs)(1)Ft,xt,x′t=0arise 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(2)Atx′t+Btxt=ftdefined on the interval I 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:(3)Atx′t+Btxt=0.In this paper, we shall consider the constant coefficient and time varying DAEs with different numerical methods. In Section2, 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 A(t) 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.
LetA∈Cn×n. If the smallest positive integer k such that (4)rankAk+1=rankAkholds it is called the index of A and is denoted by k=Ind(A). If A is nonsingular, then Ind(A)=0; else if A is singular, then Ind(A)≥1.Definition 2.
LetA∈Cn×n and Ind(A)=k. Then the matrix X∈Cn×n satisfying (1k)(5)AkXA=Ak,(2)(6)XAX=X,(3)(7)AX=XAis called the Drazin inverse of A and is denoted by X=Ad.Consider the following singular differential and corresponding difference equations(8)Atx′+Btxt=0,Atxn+1=Btxn.Supposetn+1=tn+h; h is sufficient small. Then (9)Atxtn+1-xtntn+1-tn≈Atxtn=-Btxtn.Therefore, the approximation of Ax′ is given by (10)Atxn+1=Atxn-hBtxn.Let C(t)=(A(t)-hB(t)); we can get the difference equations A(t)xn+1=C(t)xn like (8). When A(t) is a singular matrix, things can happen that are impossible when A(t)-1 exists.Definition 3.
ForA(t),B(t)∈Cn×n(t) and t0∈R, the vector c∈Cn is called a consistent initial vector associated with t0 for the equation A(t)x′(t)+B(t)x(t)=f(t) when the initial problem A(t)x′(t)+B(t)x(t)=f(t), x(t0)=c, n=1,2,…, processes at least one solution.Definition 4.
The equationA(t)x′(t)+B(t)x(t)=f(t) is said to tractable at the point t0 if the initial value problem A(t)x′(t)+B(t)x(t)=f(t), x(t0)=c has a unique solution for each consistent initial vector c associated with t0.Definition 5.
ForA(t),B(t)∈Cn×n(t) and fn∈Cn, the vector c∈Cn is called a consistent initial vector for the difference equation A(t)xn+1=B(t)xn+fn if the initial value problem A(t)xn+1=B(t)xn+fn, x0=c, n=1,2,…, has a solution for xn.Definition 6.
The difference equationA(t)xn+1=B(t)xn+fn is said to be tractable if the initial value problem A(t)xn+1=B(t)xn+fn, x0=c, n=1,2,…, has a unique solution for each consistent initial vector c.


## 3. Time Varying Differential Algebraic Equations

Lemma 7 (see [12]).
ForA(t)∈R(t)n×n, Ind(A(t))=k, and (11)gλ,t=det⁡λI-Ak+1t=g0tλn+g1tλn-1+⋯+gn-1tλ+gntis characteristic polynomial of Ak+1(t), where g0(t)=1. Then gr+1t=⋯=gn(t)=0, gr(t)≠0, and(12)Adt=-gr-1tAktAk+1tr-1+g1tAk+1tr-2+⋯+gr-2tAk+1t+gr-1tIn,where r=rank(A(t)).Next, the algorithm for the computation of the Drazin inverse of a polynomial matrix is presented as follows.Step 1.
F
i
(
t
) and gi(t) are determined by the recursive relationship as follows:(13)Fi+1t=Ak+1tFit+gitI,gi+1t=-trAk+1tFi+1ti+1,and initial conditions are(14)F1t=In,g1t=-trAk+1t.Step 2.
Fromrank(Ak(t))=rank(Ak+1(t))=r, we can see that gr+1t=⋯=gn(t)=0, gr(t)≠0. Then(15)Adt=-gr-1tAktAk+1tr-1+g1tAk+1tr-2+⋯+gr-2tAk+1t+gr-1tIn=-gr-1tAktFrt.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.
LetA(t),B(t)∈Cn×n(t). Assume that there exists λ∈C, such that (λA(t)+B(t))-1 exists, and let(16)A^λt=λAt+Bt-1At,B^λt=λAt+Bt-1Bt.Then A^λ(t)B^λ(t)=B^λ(t)A^λ(t).Proof.
If there existsλ∈C such that (λA(t)+B(t))-1 exists, then (17)λA^λt+B^λt=λλAt+Bt-1At+λAt+Bt-1Bt=λAt+Bt-1λAt+Bt=I.
Thus
B
^
λ
t
A
^
λ
t
=
I
-
λ
A
^
λ
t
A
^
λ
t
=
A
^
λ
t
-
λ
A
^
λ
t
2
=
A
^
λ
t
I
-
λ
A
^
λ
t
=
A
^
λ
t
B
^
λ
t
. Consequently, A^λ(t)B^λ(t)=B^λ(t)A^λ(t).Theorem 9.
ForA(t),B(t)∈Cn×n(t), the homogeneous difference equation(18)Atxn+1=Btxnis tractable if and only if there exists λ∈C such that (λA(t)+B(t))-1 exists.Proof.
We first prove the sufficiency. LetA^λ(t) and B^λ(t) be defined as those in (27). Clearly A(t)xn+1=B(t)xn is tractable if and only if A^λ(t)xn+1=B^λ(t)xn is tractable.
SinceλA^λ(t)+B^λ(t)=I, ∃ invertible matrix T(t)∈Cn×n(t), so that(19)T-1tA^λtTt=Ct00Nt,T-1tB^λtTt=I-λCt00I-λNt=B^1t00B^2t,where C(t) is invertible and N(t) is nilpotent of index k. Let xn=T(t)yn(t). Then the differential equation becomes(20)Ct00Ntyn+11tyn+12t=I-λCt00I-λNtyn1tyn2t,or(21)Ctyn+11t=I-λCtyn1t,Ntyn+12t=I-λNtyn2t.Since C(t) is invertible, C(t)yn+1(1)(t)=(I-λC(t))yn(1)(t) is tractable. Thus it suffices to show that the second equation of is tractable. Let k=Ind(Nt) and multiply (32) by Nk-1(t). Then (I-λN(t))Nk-1yn(2)=0, and hence Nk-1yn(2)=0. So Nk-1yn+1(2)=0. Multiply (32) again by Nk-2. Then Nk-1yn+1(2)=(I-λN(t))Nk-2yn(2), so Nk-2yn(2)=0. Continuing in this manner, we get yn(2)(t)=0 and Nyn+1(2)+(I-λN(t))yn(2)=0 is trivially tractable.
For the necessity, assume thatA(t)xn+1=B(t)xn is tractable. We need to show that there is a λ∈C such that λA(t)+B(t) is invertible. Assume that this is not true. Then λA(t)+B(t) is singular for all λ∈C. This means that, for each λ∈C, there is a vector vλ(t)∈C(t) such that (λA(t)+B(t))vλ(t)=0 and vλ(t)≠0. Let {vλ1(t),vλ2(t),…,vλs(t)} be a finite linearly dependent set of such vectors. Let xn(λi)=λinvλi(t) and let {α1(t),α2(t),…,αs(t)}⫅C(t) be such that ∑i=1s‍αivλi(t)=0, where not all αi are 0. Then zn=∑i=1s‍αixn(λi) is not identically zero and is easily seen to be a solution of (18). However, z0=∑i=1s‍λivλi(t)=0. Thus there are two different solutions, namely, zn and 0, which satisfy the initial condition x0=0. Therefore, it is not tractable at n=0, which contradicts our hypothesis. Hence, (λA(t)+B(t))-1 exists for some λ∈C.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 (λA(t)+B(t))-1 and (λA(t)-B(t))-1.Lemma 10.
Assume thatA(t),B(t)∈Cn×n(t) are such that there exists λ∈C such that (λA(t)±B(t))-1 exists. Let A^λ(t)=(λA(t)±B(t))-1A(t), B^λ(t)=(λA(t)±B(t))-1B(t). For all α,β∈C for which (αA(t)±B(t))-1 and (βA(t)±B(t))-1 exist, the following statements are true:(22)A^αdtA^αt=A^βdtA^βt,A^αdtB^αt=A^βdtB^βt,A^αtB^αdt=A^βtA^βdt,IndA^αt=IndA^βt,RA^αt=RA^βt.Theorem 11.
If the homogeneous equationA(t)xn+1=B(t)xn is tractable, then the general solution is given by(23)xn=A^tA^dtqifn=0A^dtB^tnqifn=1,2,…q∈Cnt,where A^(t)=(λA(t)-B(t))-1A(t) and B^(t)=(λA(t)-B(t))-1B(t) and λ∈C is such that (λA(t)-B(t))-1 exists. Furthermore, c∈Cn is a consistent initial vector for if and only if c∈R(A^k(t))=R(A^d(t)A(t)), where k=Ind(A^(t)). In this case this unique solution subject to x0=c is given by xn=(A^d(t)B^(t))nc, n=0,1,2,3,….Proof.
Since (3) is tractable, multiplying by (λA(t)-B(t))-1 gives the equivalent equation A^(t)xn+1=B^(t)xn. After a similarity we get, as in the proof of Theorem 9, that(24)Ct00Ntyn+11tyn+12t=I+λCt00I+λNtyn1tyn2t=B^100B^2yn1tyn2t.The difference equation is equivalent to the pair of equations(25)Ctyn+11t=B^1tyn1t,Ntyn+12t=B^2tyn2t.Since B^2 is invertible, the unique solution of the second equation of (25) was yn+1(2)(t)=B^2-k(t)NKyn+k(2)=0. But the first equation of (25) is consistent for any y0(t) and the unique solution is yn1t=C-ntB^1ty01t. In terms of the original variables we have(26)xn=Ttynt=TtC-ntB^1tn000T-1tTtI000T-1tTty01ty02t=A^dB^nq,where q=(y0(1)(t),y0(2)(t))T 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]).
LetA,B∈Cn×n. Assume that there exists λ∈C, such that (λA+B)-1 exists, and let(27)A^λ=λA+B-1A,B^λ=λA+B-1B.Then A^λB^λ=B^λA^λ.Theorem 13.
ForA,B∈Cn×n, the homogeneous difference equation(28)Axn+1=Bxnis tractable if and only if there exists λ∈C such that (λA+B)-1 exists.Proof.
The proof follows from Theorem9.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 (λA+B)-1 and (λA-B)-1.Lemma 14.
Assume thatA,B∈Cn×n are such that there exists λ∈C such that (λA±B)-1 exists. Let A^λ=(λA±B)-1A, B^λ=(λA±B)-1B. For all α,β∈C for which (γA±B)-1 and (δA±B)-1 exist, the following statements are true:(29)A^γdA^γ=A^δdA^δ,A^γdB^γ=A^δdB^δ,A^γB^γd=A^δA^δd.At last, we give the general solution of the singular difference equation.Theorem 15 (see [14]).
If the homogeneous equationAxn+1=Bxn is tractable, then the general solution is given by(30)xn=A^A^dqifn=0A^dB^nqifn=1,2,…q∈Cn,where A^=(λA-B)-1A and B^=(λA-B)-1B and λ∈C is such that (λA-B)-1 exists. Furthermore, c∈Cn is a consistent initial vector if and only if c∈R(A^k)=R(A^d(t)A), where k=Ind(A^). In this case this unique solution subject to x0=c is given by xn=(A^dB^)nc,n=0,1,2,3,….Proof.
The proof follows from Theorem11


### 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.Thes-stage implicit Runge-Kutta method applied to the general nonlinear DAE of the form (2) is defined by(31)0=Fti,Xi,Ki,(32)ti=tn-1+cih,i=1,2,…,s,(33)Xi=xn-1+h∑j=1saijKj,(34)xn=xn-1+h∑i=1sbiKi.The method is often denoted by the shorthand notation or Butcher diagram: (35)cAb^T=c1a11a11⋯a1sc2a21a22⋯a2s⋮⋮⋮⋱⋮csas1as2⋯assb1b2⋯bs.


### 4.2. Padé Approximation [15]

A differential algebraic equation has the form(36)Ft,x,x˙=0with initial values (37)xt0=x0,x˙t0=x1,where F and x are vector functions for which we assumed sufficient differentiability.We assume the solution has the form(38)x=x0+x1t+et2,where e is a vector function which is the same size as x0 and x0′. Substituting (38) into (36) and neglecting higher order term, we have the linear equation of e in the form(39)Ae=B,where A and B are constant matrixes. Solving (39), the coefficients of t2 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].


## 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.Thes-stage implicit Runge-Kutta method applied to the general nonlinear DAE of the form (2) is defined by(31)0=Fti,Xi,Ki,(32)ti=tn-1+cih,i=1,2,…,s,(33)Xi=xn-1+h∑j=1saijKj,(34)xn=xn-1+h∑i=1sbiKi.The method is often denoted by the shorthand notation or Butcher diagram: (35)cAb^T=c1a11a11⋯a1sc2a21a22⋯a2s⋮⋮⋮⋱⋮csas1as2⋯assb1b2⋯bs.


## 4.2. Padé Approximation [15]

A differential algebraic equation has the form(36)Ft,x,x˙=0with initial values (37)xt0=x0,x˙t0=x1,where F and x are vector functions for which we assumed sufficient differentiability.We assume the solution has the form(38)x=x0+x1t+et2,where e is a vector function which is the same size as x0 and x0′. Substituting (38) into (36) and neglecting higher order term, we have the linear equation of e in the form(39)Ae=B,where A and B are constant matrixes. Solving (39), the coefficients of t2 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:(40)Ax˙+Bx=0,where(41)A=23210-2000,B=181410012-27-21-15,and we take h=0.1; the analytical solution of (40) is x(t)=(1/6+1/6e2t/3,-1/3+4/3e2t/3,1/6-13/6e2t/3)T.(1) Drazin Inverse. Suppose tn+1=tn+h; h is sufficiently small. Then (42)Axtn+1-xtntn+1-tn≈Axtn=-Bxtnso on the point tn+1=tn+h the approximate value of Ax′ is given by (43)Axn+1=Axn-hBxn.Now we use the difference rule to solve the equation. Set(44)C=A-hB=0.20001.60001.00001.0000-0.1000-2.20002.70002.10001.5000,and then (45)A+C=2.20004.60003.00002.0000-0.1000-4.20002.70002.10001.5000.det⁡A+C=-33.48≠0, so λ=1 such that (λA+C)-1 exists; we can see the homogeneous equation Axn+1=Cxn is tractable. For -A-C=-(A+C), we have λ=-1 such that (λA-C)-1 exists. Thus (46)A^=-A-CA=0.50000.77690.5538-1.0000-1.2849-0.56990.50000.4005-0.1989,C^=-A-CC=-1.5000-0.7769-0.55381.00000.28490.5699-0.5000-0.4005-0.8011.The eigenvalues of A^d are 0, -0.5, and -0.4839, so that A^d could be computed by Theorem  7.5.2 in [13]: (47)A^d=2.00003.11482.2296-4.0000-5.0815-2.16302.00001.5074-0.9852.According to Theorem 11, the general solution is given by (48)xn=-1.0000-1.5556-1.11112.00002.55561.1111-1.0000-0.77780.4444qifn=0-1.0000-1.5593-1.11852.0002.52591.0519-1.000-0.72960.5407nqifn=1,2,…,where q=x(0)=(1/3,1,-2)T, which is a consistent initial vector.(2) Radau IIA Methods. We consider Radau IIA as follows: (49)13512-112134143414The 2-stage implicit Runge-Kutta method is applied to (40):(50)AKn,1+Bxn+h512Kn,1-112Kn,2=0,AKn,2+Bxn+h34Kn,1-14Kn,2=0.We can rewrite (50) as follows:(51)A+512hB-112hB34hBA+14hBKn,1Kn,2=-Bxn-Bxn,where the coefficient matrix is nonsingular. So, we can compute the numerical solutions by the following equation: (52)xn+1=xn+h34Kn,1+14Kn,2.(3) Padé Approximation. From initial values x(0)=(1/3,1,-2)T, the solutions of (40) can be supposed as(53)x1t=13+e1t,x2t=1+e2t,x3t=-2+e3t.Substituting (53) into (40) and neglecting higher order terms, we have(54)2e1+3e2+2e3+Ot=0,e1-2e3+Ot=0,9e1+7e2+5e3t=0.The linear equation can be given in the following:(55)23210-2975e1e2e3=000.Solving this linear equation we have e=(1/9,8/9,-13/9)T; then the solutions of (40) can be supposed as(56)x1t=13+19t+e1t,x2t=1+89t+e2t,x3t=-2-139t+e3t.In the same manner, substituting (56) into (40) and neglecting higher order terms, then we have(57)4e1+6e2+4e3t+Ot2=0,2e1-4e3-2t+Ot2=0,-27e1-21e2-15e3t2=0.Similar to (53), solving this linear equation, we have e=(1/27,8/27,-13/27)T; then the solutions of (40) can be supposed as(58)x1t=13+19t+127t2+e1t3,x2t=1+89t+827t2+e2t3,x3t=-2-139t-1327t2+e3t3.Repeating the above procedure we have (59)x1t=13+19t+127t2+2243t3+1729t4+210935t5+298415t6+42066715t7+16200145t8+2167403915t9+Ot10,x2t=1+89t+827t2+16243t3+8729t4+1610935t5+1698415t6+322066715t7+86200145t8+16167403915t9+Ot10,x3t=-2-139t-1327t2-26243t3-13729t4-2610935t5-2698415t6-522066715t7-136200145t8-26167403915t9+Ot10.The power series can be transformed into the following Padé series: (60)v1t=54=1/3+1/81t+4/243t2+1/1701t3+1/15309t4+1/688905t51-8/27t+1/27t2-4/1701t3+1/15309t4,v2t=54=1+16/27t+17/243t2+44/5103t3+19/45927t4+8/688905t51-8/27t+1/27t2-4/1701t3+1/15309t4,v3t=54=-2-23/27t-31/243t2-67/5103t3-32/45927t4-13/688905t51-8/27t+1/27t2-4/1701t3+1/15309t4.If the power seriesx1(t),x2(t),x3(t) converges very fast then x1(t),x2(t),x3(t) can be eliminated in Padé series.In Tables1–5, exact solutions, numerical solutions, and errors are illustrated.Table 1
Exact solution and numerical solution ofx(1)(n).
t (h = 0.1)
x
(
1
)
(
n
)
x
n
(
1
) of Drazin inv
x
n
(
1
) of Radau
x
n
(
1
) 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.4912890064Table 2
Exact solution and numerical solution ofx(2)(n).
t (h = 0.1)
x
(
2
)
(
n
)
x
n
(
2
) of Drazin inv
x
n
(
2
) of Radau
x
n
(
2
) 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.2636453884Table 3
Exact solution and numerical solution ofx(3)(n).
t (h = 0.1)
x(3)(n)
x
n
(
3
) of Drazin inv
x
n
(
3
) of Radau
x
n
(
3
) 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.0534237552Table 4
ErrorsEn≔xn-xn2.
t (h = 0.1)
E
n of Drazin inverse
E
n of Radau
E
n of Padé approximation
0.1
5.9575917383
e
-
003
7.6143396932
e
-
007
5.9160779268
e
-
010
0.2
1.2916314922
e
-
002
1.6276920132
e
-
006
4.5825760741
e
-
010
0.3
2.0777743067
e
-
002
2.6097631713
e
-
006
2.4494899455
e
-
010
0.4
2.9630445273
e
-
002
3.7196628261
e
-
006
5.3851633972
e
-
010
0.5
3.9572170567
e
-
002
4.9700993038
e
-
006
4.6904157928
e
-
010
0.6
5.0709583610
e
-
002
6.3750949076
e
-
006
6.5574377970
e
-
010
0.7
6.3158869695
e
-
002
7.9504631349
e
-
006
8.2462103179
e
-
010
0.8
7.7046497345
e
-
002
9.7126908162
e
-
006
7.5498350599
e
-
010
0.9
8.3510065024
e
-
002
1.1679886007
e
-
005
5.9160780206
e
-
010
1.0
9.1969922436
e
-
002
1.3872451324
e
-
005
6.4031244204
e
-
010Table 5
Errors of Drazin inverse and Radau IIA.
t (h = 0.01)
E
n of Drazin inverse
E
n of Radau IIA
0.01
1.1228248162
e
-
004
1.0088260628
e
-
003
0.02
6.3035622437
e
-
004
8.6145976589
e
-
003
0.03
5.5249099531
e
-
004
1.2036909993
e
-
003
0.04
3.9265521061
e
-
004
1.3035988488
e
-
003
0.05
1.7558562753
e
-
004
5.5498031923
e
-
003
0.06
1.7313059801
e
-
004
6.7699688079
e
-
003
0.07
1.2203796056
e
-
004
7.1281533560
e
-
003
0.08
5.0551368477
e
-
004
9.6806736561
e
-
003
0.09
6.6468312281
e
-
004
9.9061602809
e
-
003
0.1
6.1145689285
e
-
004
1.3892356616
e
-
003Example 2 (see [16]).
We consider the following differential algebraic equation:(61)Atx′t+Btxt=0,where (62)A=1000010000100000,B=10-t1-11-t2t-t3t2-10t-1t-1,and we take h=0.01; the analytical solution of (40) is x(t)=(e-t,te-t,e-t,tet)T.(1) Drazin Inverse. The classical four-order Runge-Kutta method is applied to (40); we can get (63)Atxn+1=Atxn+hk1+2k2+2k3+k46,where (64)k1=-Btxn,k2=-Bt+hk12,k3=-Bt+hk22,k4=-Bt+hk3.Thus (65)Atxn+1=At-h+h22+h36+h424Btxn.Now we use the difference rule to solve the system of (40). Set (66)Ct=At-h+h22+h36+h424Bt=0.894800.1052t-0.10520.10520.89480.1052t2-0.1052t0.1052t3-0.1052t21.10520-0.1052t0.1052-0.1052t0.1052;set λ=-1; we get det(-A(t)-C(t))=0.7925≠0, so we have λ=-1 such that (λA(t)-C(t))-1 exists. Thus (67)A^t=-At-CtA=-0.5294-0.0294t0.0294000.0294-0.0294t2-0.5294+0.0294t000.0015t2+0.0264t3-0.0264t2-0.4750-0.0294-0.5294t+0.0015t3+0.0264t40.5294-0.0264t3-0.475t0,C^t=-At-CtC=-0.4975+0.0025t-0.002500-0.0025+0.0025t2-0.4975-0.0025t00-0.0025t30.0025t2-0.502500.0025+0.5025t-0.0025t4-0.5025+0.0025t30.4975t-1so that A^d could be computed by finite algorithm (13)–(15): (68)A^d=-1.8948+0.1052t-0.105200-0.1052+0.1052t2-1.8948-0.1052t00-0.1052t30.1052t2-2.105200.1052-1.8948t-0.1052t41.8948+0.1052t3-2.1052t0.According to Theorem 15, the general solution is given by (69)xn=1-2.8544t2.8544002.8544-2.8544t21+2.8544t001.426t2+2.5692t3-2.5692t210-2.8544+t+1.426t3+2.5692t4-1-2.5692t3t0qifn=00.8948-0.1052t0.1052000.1052-0.1052t20.1052t00-1.8334t+0.1052t3-0.1052t21.10520-0.1052+0.8948t-1.8334t3+0.1052t4-0.8948-0.1052t31.1052t0qifn=1,2,…,where q∈Cn(t).Letq∈R(A^(t))k; we can compute the numerical solution of (40) by (33).(2) Radau IIA. We consider the Radau IIA as follows: (70)13512-112134143414Radau IIA is applied to (40): (71)AKn,1+Bxn+h512Kn,1-112Kn,2=0,AKn,2+Bxn+h34Kn,1-14Kn,2=0.We can rewrite (50) as follows:(72)A+512hB-112hB34hBA+14hBKn,1Kn,2=-Bxn-Bxn,where the coefficient matrix is nonsingular. So, we can compute the numerical solutions by the following equation: (73)xn+1=xn+h34Kn,1+14Kn,2.See Table5 for errors of the numerical results. Error En is defined by En≔xn-xn2.


## 6. Conclusion

The fundamental goal of this study has been to construct approximations to numerical solutions of linear DAEs. From Table4, 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.



---

*Source: 1871590-2017-06-08.xml*