Abstract
This paper presents a generalized high order block hybrid -step backward differentiation formula (HBDF) for solving stiff systems, including large systems resulting from the semidiscretization parabolic partial differential equations (PDEs). A block scheme in which two off-grid points are specified by the zeros of the second degree Chebyshev polynomial of the first kind is examined for convergence, and stabilities. Numerical simulations that illustrate the accuracy of a Chebyshev based method are given for selected stiff systems and partial differential equations.
1. Introduction
We consider the first order differential equation where , , , satisfies a Lipschitz condition, and the eigenvalues of the Jacobian have negative real parts (see [1]). It is well known that the system (1) is better handled by methods with larger stability intervals. In particular, -stable methods are of great importance. However, for very large systems arising from the semidiscretization of parabolic PDEs, -stable methods converge very slowly to the exact solution. Hence, we seek methods which are at least -stable for efficiently solving (1) when the system is very large (see Cash [2]).
The development of continuous methods has been the subject of growing interest due to the fact that continuous methods enjoy certain advantages, such as the potential for them to provide defect control (see Enright [3]) as well as having the ability to generate additional methods, which are combined and applied in block form (see Onumanyi et al. [4, 5], Akinfenwa et al. [6], and Jator [7]).
The majority of block methods which are due to Shampine and Watts [8], Chartier [9], Rosser [10], and Chu and Hamilton [11] are generally implemented in the predictor-corrector mode. In this paper, we adopt a different approach where the solution is simultaneously provided in each block (see Jator et al. [12], Jator [7, 13]) without the use of predictors from other methods. It is unnecessary to make a function evaluation at the initial part of the new block since at all blocks (with the exception the first block) the first function evaluation is already available from the previous block.
The paper is structured as follows. Section 2 states the generalized block HBDF, its derivation from continuous approximation, and how to generate specific members of the scheme. The HBDFs are bundle as main and additional methods, a concept that is due to Brugnano and Trigiante [14]. In Section 3 we study the stability of the schemes with emphasis on a Chebyshev based member. In Section 4 a numerical algorithm for the block method is provided and implemented on selected test problems. Finally, the conclusion of the paper is discussed in Section 5.
2. Hybrid Backward Differentiation Type Formulas
2.1. Step HBDF
We define a -step HBDF (main method) for the numerical solution of (1) on the interval as where is the step number, , , , , , are constant coefficients, and are off-grid points which cannot be integers. We note that , , and . Conventionally the method (2) is applied in the predictor-corrector mode which involves the use of starting values and predictors from other methods. In this paper, a different approach that involves a block-by-block implementation is adopted, whereby additional methods are generated from a continuous scheme and combined with (2) to form a block method which is then applied to (1) without the need for starting values and predictors. In what follows, we state the block -step HBDF for solving (1).
2.2. Block HBDF
We define the block -step HBDF as where the first two members of (3) are additional methods provided by the continuous scheme constructed in the next subsection; , , , , , , , , , , , are constant coefficients, which are determined via the procedure given in the next subsection. The numerical solution is an approximation to the analytical solution , , , and . In order to facilitate the analysis of (3) we express it in block form as follows: where where and , and , , , are by matrices whose entries are given by the coefficients of (3). We note that is the number of blocks. The derivation given in the next subsection will be used to produce (4).
2.3. Derivation of Continuous HBDF
We provide here a generalized framework on which the schemes stated in the preceding section are produced. Let be the linear space of polynomials of degree and let be a given polynomial. On the interval from to , the exact solution for (1) is approximated by where and are continuous coefficients.
Let the following conditions be satisfied: Further, suppose that (6) is defined by the polynomial basis functions where the constants , , and are uniquely determined.
Substituting (8) into (6) we get which simplifies to and can be expressed as where
By imposing the conditions (7) on (11), we obtain a system of equations, which can be expressed as where , , and is a matrix of dimension defined as where , , are basis functions. The elements of the vector are obtained using Cramer’s Rule as where is obtained by replacing the th column of by . The values of are substituted into (11) to obtain the continuous approximation which is used to generate the main and additional methods.
Remark 1. The coefficients of (6) are given by the explicit representation (16), since the two methods are equivalent. We note that the explicit representation given by (16) is only used for derivation purposes and is not recommended to be used numerically since its implementation will be expensive.
2.4. Generating Specific Members of the Method (3)
We illustrate here how to generate specific members of the method (3) and in particular present the schemes for the case where the off-grid points are chosen as the zeros of Chebyshev’s polynomial. If we choose and for notational simplicity let , then on the interval from to , the approximation (6) becomes with first derivative given by where , , and are continuous coefficients uniquely determined as and where , , , , and . We also introduce the scale variable for convenience only by letting .
Suppose that the off-grid points are specified by the zeros of the second degree Chebyshev polynomial of the first kind given by and ; then (18) yields three additional methods when evaluated at which are combined with the main method obtained by evaluating (17) at to get a Chebyshev based method of the form (3) given as whose stability properties are discussed in the next section.
3. Stability of Block HBDF
In this section, we discuss the stability of the block HBDF (4) with emphasis on the Chebyshev based method (20). We pause to remark that the method (20) is a specific member of the class of methods (3) and therefore a similar analysis can be generated for other members.
Suppose the local truncation error of the block method (4) is defined as where We recall here that the block method (4) is characterized by the matrices , , , whose entries are given by the coefficients of (3).
Definition 2. The block method (4) has order of consistency .
Thus, we expand the terms in (21) as a Taylor series about the point and also assume that the constant coefficients for , leads to the local truncation error satisfying the maximum norm .
Theorem 3. The block method (4) is zero-stable.
Proof. Zero-stability is concerned with the stability of a difference system in the limit as tends to zero. Thus, as and noting that the first columns of are all zeros with nonzero elements occurring only in the last column, the normalized first characteristic polynomial of the method (4) is given by is a constant and is the identity matrix. Following Fatunla [15], the block method (4) is zero-stable, since from (23), satisfy , , and for those roots with , the multiplicity does not exceed 1.
In particular, the Chebyshev based method (20) is zero-stable as a consequence of Theorem 3, since our analysis shows that its associated characteristic polynomial is given by . Thus, the root condition is satisfied, since for , , , and for those roots with , the multiplicity does not exceed 1.
Lemma 4. The block method (4) is convergent.
Proof. Following Henrici [16] that convergence = consistency + zero-stability, the proof is a consequence of Definition 2 and Theorem 3.
The stability of the method is studied next. The application of (4) to the test equation gives the following discrete problem: where is the stability matrix.
Definition 5. The block method (4) is said to be (i) -stable if for all , has a dominant eigenvalue such that ; moreover, since is a rational function, the real part of the zeros of must be negative, while the real part of the poles of must be positive; (ii) -stable if for all , has a dominant eigenvalue such that ; (iii) -stable if it is -stable and as ; and (iv) -stable if it is -stable and as .
Definition 6. Let the stability region of (4) be ; then the method (4) is said to be -stable, , if .
Theorem 7. The block HBDF (4) is -stable and -stable.
Proof. The eigenvalues of the stability matrix in (24) are rational functions of . It can be easily shown that the dominant eigenvalue has the form , where the degree of the polynomial is greater than that of and . Applying Definitions 5 and 6 with completes the proof.
Theorem 8. The Chebyshev based method (20) is -stable and -stable.
Proof. The dominant eigenvalue for the Chebyshev based method is given by The proof follows from Theorem 7.
Moreover, the stability region for the Chebyshev based method is given in Figure 1 showing the zeros and poles of the dominant eigenvalue .
4. Numerical Application
4.1. Block Algorithm
Given the interval where for , the method (3) is implemented in block form (4) as follows.
Step 1. Let ; the method (3) reduces to the Chebyshev based method (20), and for , and , the values of ( is the transpose) are simultaneously obtained over the subinterval since is known from the system of ODEs.
Step 2. For , , the values of are simultaneously obtained over the subinterval since the values of are known from the previous block.
Step 3. For , the process is continued to generate the vectors for the approximate solutions of the ODE on subintervals .
4.2. Numerical Examples for Stiff ODEs
We will implement the block algorithm described above using the Chebyshev based method (20) for three different stiff systems given in the literature. The first example shows the performance of the block method on a linear stiff system. The next two examples are chosen to demonstrate the performance of the block method on a and a nonlinear stiff system.
Example 1. Consider the following linear system on the range : , , and with initial conditions , , and , respectively. The exact solution of the system is given by , , and .
This problem has also been solved by Amodio and Mazzia [17] using the boundary value methods (BVMs), implicit Adams methods (Ad-IVMs), and backward differentiation formulas (BDFs). The results for their order 4 methods are reproduced in Table 1 and compared with the results given by the HBDF which is also of order 4. It is seen from Table 1 that the HBDF performs better than those in [17]. In all cases, the calculated rate of convergence of each method is consistent with the theoretical order () behavior of each method. In particular, it is obvious from Table 1 that HBDF exhibits an order 4 behavior, since on halving the step size, the error is reduced by a factor of about 18, which can be expressed approximately as . Thus, for this example, HBDF is superior in terms of accuracy. We note that the maximum relative errors for the first component (Errs) displayed in Table 1 are computed as max .
Example 2. We consider the following IVP which was solved by Vaquero and Vigo-Aguiar, [18]: and with initial conditions and , respectively. The exact solution is ; .
The results for this problem obtained by Vaquero and Vigo-Aguiar [18] using an exponentially fitted Gauss (EF-Gauss-2s) and Gauss-2s methods of order 4 are reproduced in Table 2 and compared with the results given by HBDF which is also of order 4. It is seen from Table 2 that the HBDF performs better than those in [18].
Example 3. We consider the given nonlinear system on the range : True solution is (see Jeltsch [19]).
This example was chosen to show the performance of the HBDF on a nonlinear stiff system. In Table 3, the accuracy of the HBDF is measured by the maximum computing absolute errors at .
4.3. Numerical Examples for Parabolic PDEs
In this section we test the block algorithm with Chebyshev based scheme (20) on large systems derived from parabolic pdes. Consider the homogeneous partial differential equation subject to the initial/boundary conditions , , , and .
Suppose that (28) is discretized using the method of lines in such as way that the resulting system of ordinary differential equations is stable (see Lambert [20], Ramos and Vigo-Aguiar [21], and Cash [2]). Let be the exact solution of the resulting semidiscrete problem where , , , and The problem (29) is now a system of ordinary differential equations on which the block algorithm can be applied. We remark that if is small, the system (29) is stiff.
Example 4. Consider the given PDE (see Cash [2]) The exact solution .
The output in Table 4 compares the performance of the block method with those given in Cash [2]. It is clear from the table that the block method is the more accurate.
Example 5. Consider the given stiff parabolic equation (see Cash [2]) The exact solution .
Cash [2] notes that as increases, equations of the type given in Example 5 exhibit characteristics similar to model stiff equations. Hence, methods such as the Crank-Nicolson method which are not -stable are expected to perform poorly. The block HBDF is -stable and hence performs well on this problem. Table 5 displays the results for and a range of values for .
5. Conclusion
In this paper, we have proposed a generalized high order block hybrid -step HBDF for the numerical solution of stiff systems, which includes large systems arising from semidiscretized parabolic PDEs. A specific scheme based on the zeros of the second degree Chebyshev polynomial of the first kind was generated and its stability established. Applications of the block algorithm to stiff ODEs and PDEs demonstrate the accuracy of the method. The block algorithm is facilitated by the availability of the continuous method which has the ability to generate additional methods, which are combined and applied in block form. We remark that the block-by-block approach has stability and accuracy advantages over the step-by-step methods which are generally applied in a predictor-corrector mode. Particularly, the block-by-block approach gives smaller global errors at the end of the interval compared to those produced by the step-by-step methods due to the fact that the accumulation of error at each step is inherent in the step-by-step methods. Our future research will be focused on applying continuous methods to delay and differential algebraic equations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.