BLOCK METHOD FOR NUMERICAL SOLUTION OF DIFFERENTIAL-ALGEBRAIC EQUATIONS WITH HESSENBERG INDEX-3

A Block Second Derivative Formula (BSDF) of order seven for the numerical solution of Hessenberg Differential Algebraic Equations (DAEs) of index 3 is presented. This is achieved by constructing a Continuous Second Derivative Formula (CSDF) used to obtain the main and additional methods which are combined to form a single block of methods that simultaneously provide the approximate solutions for the Hessenberg DAEs of index 3. The error constant and stability properties of the BSDF are discussed. The performance of the method is demonstrated on some numerical examples to show the accuracy and efficiency of the method. AMS Subject Classification: 65L05, 65L06, 65L08


Introduction
Many important mathematical problems can be modeled as systems of Differential-Algebraic Equations (DAEs).These problems have a wide range of applications in various branches of science and engineering.These include mechanical or multibody systems, chemical processes, optimal control, electric circuit design, analytical surveys and dynamical systems.In the literature, some numerical methods have been developed for the solution of DAEs such as the BDF (see [2], [3], [9]), implicit Runge-Kutta methods [3], Pade and Chebysev approximation methods (see [4] [5] [6]) and variational iterative method [12].These methods are only directly suitable for low-index problems and often require that the problem, have special structure.Although many important applications can be solved by these methods, there is a need for more efficient, accurate and reliable techniques.
A system of ordinary differential equations (ODEs) with algebraic constraints which can be written in the general form as is called differential algebraic equation, where ∂F ∂y ′ may be singular and the structure and rank of the of this Jacobian matrix may depend on the solution y(t).An important special case of ( 1) is the semi-explicit DAE of the form whose index is 1 if ∂g ∂z is nonsingular, since in principle differentiating f 2 (y(t), z(t)) = 0, once yields z ′ .Thus the distinction between the differential variables y(t) and algebraic variables z(t).The algebraic variables may be less smooth than the differential variables by one derivative.In the general case, each component of y may contain a mix of differential and algebraic components, which makes the numerical solutions of such high-index problems harder and with higher risk.
Most of the higher-index problems encountered in practice can be expressed as a combination of more restrictive structures of ODEs coupled with constraints.In such systems the algebraic and differential variables are explicitly identified and all the algebraic variables may be eliminated using the same number of differentiations.These type of DAEs are called Hessenberg DAEs.
Definition 1.1 (Hessenberg Index-1) Here the Jacobian matrix function f 2z is assumed to be nonsingular for all t.This is also often referred to as a semi-explicit index-1 system.Definition 1.2 (Hessenberg Index-2) the variables z of the algebraic part is absent from the constraint and the product of Jacobians f 2y f 1z is nonsingular for all t.Thus, all algebraic variables play the role of index-2 therefore this is a pure index-2 DAE.Definition 1.3 (Hessenberg Index-3) In this case the product of three matrix functions f 3w f 2y f 1z is nonsingular.The index of a Hessenberg DAE is found, as in the general case, by differentiation.However, only the algebraic constraints must be differentiated.
The aim in this paper is to develop a Block Second Derivative formula (BSDF) of order k+2 for the solution of Hessenberg DAE of the form (5) that will be efficient, reliable and accurate.Block methods were first introduced by Milne [14] for use only as a means of obtaining starting values for predictorcorrector algorithms and has since then been developed by several researchers (see [ [8], [16], [17], [18]).Recently, Akinfenwa and Okunuga [1] and Naghmeh et.al.[15] developed block method for the solution of semi explicit index 1 DAEs This paper as in [1] presents a block method which preserves the Runge-Kutta traditional advantage of being self-starting and efficient, since it requires m function evaluations per integration step, where m is the number of functions in the block method.

Construction of BSDF
In this section, we construct the continuous scheme of the main and additional methods known as the second derivative formula (SDF) and are combined to form the BSDF on the interval from t n to t n+k = t n + kh where h is the chosen step-length and k is the step number of the method.We assume that the exact solution y(t) on the interval [t n , t n+k ] is locally represented by Y (t) given by l j are unknown coefficients to be determined, and ϕ j (t) are polynomial basis function of degree s + r − 1 such that the number of interpolation points s and the number of distinct collocation points r are respectively chosen to satisfy 1 ≤ s < k, k = 5, and r > 0. The proposed class of methods is thus constructed by imposing the following conditions k+2 j=0 l j t j n+i =y n+i , i = 0, (7) assuming that y n+i = Y (t n + ih), denote the numerical approximation to the exact solution y(t n+i ) f n+i = Y ′ (t n + ih), y n+j ),denote the approximation to y ′ (t n+i ) n is the grid index.It should be noted that equation ( 7),( 8) and ( 9) lead to a system of k + 3 equations which must be solved to obtain the coefficients l j , j = 0, . . ., 7. The continuous second derivative formula (CSDF) is then obtained by substituting the values of l j into equation ( 6).After some algebraic reduction, the CSDF yields the expression in the form ( 10) where α j (t), j = 0, β j (t), j = 0, 1, . . ., 5 and γ 5 (t) are continuous coefficients.
The CSDF (10) is then used to generate the main and additional discrete formulas by evaluating at points t = t n+5 and t = t n+1 , . . ., t n+4 ].These discrete formulas given in ( 11) and ( 12) are then combined and implemented simultaneously as a single block formula for the solution of the Hessenberg index 3 DAEs.

Analysis of the BSDF
The SDFs can be represented by a matrix finite difference equation in block form as where: for ω = 1, . . .π N , π N is number of partition on [a,b] and n = 0, 5, . . ., N − 5.
And A (1) is an identity matrix of dimension 5, while the matrices A (0) , B (1) , B (0) and C (1) are 5 by 5 matrices whose entries are given by the combined coefficients of the methods obtained in (10) given as follows: Following Fatunla [7] and Lambert [13], the local truncation error associated with each of the method in the BSDF can be defined by a linear difference operator Assuming that y(t) is sufficiently differentiable, we can write the terms in ( 14) as a Taylor series about the point t to obtain the expression where the constant coefficients C p , p = 0, 1, 2, . .., are given as follows: . . .
According to Henrici [11], we say that the methods in (11) have a maximal order of accuracy p if Therefore, C p+1 is the error constant and C p+1 h p+1 y (p+1) (t n ) the principal local truncation error at the point t n .Hence from our calculations the order of the BSDF is seven and the local truncation error is given as ) where T is transpose.

Zero Stability
The zero stability of the method is concerned with the stability of the difference system in the limit as h tends to zero (see [7]).Thus, as h −→ 0 the difference system (13) tends to Whose first characteristics polynomial ρ(R) given by .

Consistency and Convergence
We note that the BSDF is consistent as it has order p > 1 Since the BSDF ( 13) is zero stable.According to Henrici [11], Convergence = zero stability + consistency.Hence the BSDF converges.

Linear stability
The linear stability properties of the BSDF is discussed in the spirit of Hairer and Wanner [10] and determined by expressing it in the form (13) and applying the test problem where the matrix Q(z) is given by The matrix Q(z) has eigenvalues {ξ 1 , ξ 2 , . . ., ξ 5 } = {0, 0, . . ., ξ 5 }, where the dominant eigenvalue ξ 5 is the stability function ξ(z) : C → C which is a rational function with real coefficients given by 6(420 + 900z + 850.z 2 + 450z 3 + 137z 4 + 20z 5 ) 2520 − 7200.z+ 9600z 2 − 7800z 3 + 4197z 4 − 1490z 5 + 300z 6 The stability domain of the method (or region of absolute stability) S is defined as Specifically, when the left-half complex plane is contained in S, the method is said to be A-stable.In Figure 1, the plot with the rectangles representing the zeros and plus signs representing the poles of (19).The plot in white represents the stability region which corresponds to the stability function (19).Clearly, from Figure 1, it is obvious that our method is A-stable since according to Hairer and Wanner [10] there is no pole of the stability function (19) in the left half complex plane.
For n = 1, ω = 2, (y 6 , y 7 , . . ., y 10 ) T are simultaneously obtained over the sub-interval [t 5 , t 10 ], as y 5 is known from the previous block.Hence, the subintervals do not over-lap.It should be noted that for linear problems, the code uses Gaussian elimination while the Newton's method is used for nonlinear problems.

Numerical Examples
This section deals with some numerical experiments, executed in our written code in MATHEMATICA 8 to illustrate the result derived in the previous sections.
Example 5.1.: Consider the differential-algebraic equations with Hessenberg index-3 given by The exact solution is given as: Example 5.2.: Our second example is the differential-algebraic equations with Hessenberg index-3 given by whose exact solution is given as: The results for the two problems are tabulated for h = 0.1 and compared with the results in [12]

Conclusion
We have proposed a block second derivative formula of order seven for the solution of differential-algebraic equations DAEs with Hessenberg index-3.The algorithms are self-starting, consistent and provide good accuracy.Numerical examples using the BSDF showed that the method is accurate and efficient for the class of problems being considered.

Table 2 :
comparison of methods for y *

Table 4 :
comparison of methods for y *

Table 6 :
comparison of methods for y *