# Block Hybridk-Step Backward Differentiation Formulas for Large Stiff Systems

Authors: S. N. Jator; E. Agyingi
Journal: International Journal of Computational Mathematics
 (2014)

---

## Abstract

This paper presents a generalized high order block hybridk-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, L and A stabilities. Numerical simulations that illustrate the accuracy of a Chebyshev based method are given for selected stiff systems and partial differential equations.


---

## Body

## 1. Introduction

We consider the first order differential equation(1)y′(t)=f(t,y(t)),y(t0)=y0,t∈[t0,tN],
where N∈Z+, f:R×Rm→Rm, y,y0∈Rm, f satisfies a Lipschitz condition, and the eigenvalues of the Jacobian (∂f/∂y) 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, A-stable methods are of great importance. However, for very large systems arising from the semidiscretization of parabolic PDEs, A-stable methods converge very slowly to the exact solution. Hence, we seek methods which are at least L0-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. Section2 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.k-Step HBDF

We define ak-step HBDF (main method) for the numerical solution of (1) on the interval [tn,tn+k] as
(2)∑j=0k‍αjyn+j+∑j=1k‍ανjyn+νj=hβkfn+k,
where k∈Z+ is the step number, tn=t0+nh, h=(tN-t0)/N, n=0,…,N, αj, ανj, βk are constant coefficients, and νj∈(0,k) are off-grid points which cannot be integers. We note that yn+j≈y(tn+j), yn+νj≈y(tn+νj), and fn+k=f(tn+k,yn+k). 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 k-step HBDF for solving (1).


### 2.2. Block HBDF

We define the blockk-step HBDF as
(3)∑j=0k-1‍α¯jyn+j+∑j=1k‍α¯νjyn+νj=hβ¯kfn+k+hβ¯ifn+i,i=1,…,k-1,∑j=0k-1‍α^jyn+j+∑j=1k‍α^νjyn+νj=hβ^kfn+k+hβ^νifn+νi,i=1,…,k,∑j=0k‍αjyn+j+∑j=1k‍ανjyn+νj=hβkfn+k,
where the first two members of (3) are additional methods provided by the continuous scheme constructed in the next subsection; αj, ανj, βνi, βk, α¯j, α¯νj, β¯νi, β¯k, α^j, α^νj, β^νi, β^k are constant coefficients, which are determined via the procedure given in the next subsection. The numerical solution yn+j is an approximation to the analytical solution y(tn+j), yn+νj≈y(tn+νj), fn+νj=f(tn+νj,yn+νj), and fn+j=f(tn+j,yn+j). In order to facilitate the analysis of (3) we express it in block form as follows:
(4)A1Wμ=A0Wμ-1+h[B1Fμ+B0Fμ-1],
where
(5)Wμ=(yn+ν1,…,yn+νk,yn+1,…,yn+k)T,Wμ-1=(yn-ν1,…,yn-νk,yn-k+1,…,yn)T,Fμ=(fn+ν1,…,fn+νk,fn+1,…,fn+k)T,Fμ-1=(fn-ν1,…,fn-νk,fn-k+1,…,fn)T,
where μ=1,…,Ϝ and n=0, k,…,N-k and Ai, Bi, i=0,1, are 2k by 2k 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ΥJ be the linear space of polynomials of degree ≤J and let U(t)∈ΥJ be a given polynomial. On the interval from tn to tn+k=tn+kh, the exact solution y(t) for (1) is approximated by
(6)U(t)=-∑j=0k-1‍αj(t)yn+j-∑j=1k‍ανj(t)yn+νj+hβk(t)fn+k,
where αvj(t) and βk(t) are continuous coefficients.Let the following conditions be satisfied:(7)U(tn+vj)=yn+vj,j=1,2,…,k,U(tn+j)=yn+j,j=0,1,2,…,k-1,U′(tk)=fn+k.
Further, suppose that (6) is defined by the polynomial basis functions
(8)αj(t)=∑i=02k‍αi+1,jPi(t),j={0,1,…,k-1},αvj(t)=∑i=02k‍αi+1,jPi(t),j={1,…,k},hβk(t)=∑i=02k‍hβi+1,kPi(t),
where the constants αi+1,j, αi+1,vj, and hβi+1,k are uniquely determined.Substituting (8) into (6) we get
(9)U(t)=-∑j=0k-1‍∑i=02k‍αi+1,jPi(t)yn+j-∑j=1k‍∑i=02k‍αi+1,jPi(t)yn+νj+h∑i=02k‍hβi+1,kPi(t)fn+k,
which simplifies to
(10)U(t)=∑i=02k‍{-∑j=0k-1‍αi+1,jyn+jmm-∑j=1k‍αi+1,jyn+νj+hβi+1,kfn+k}Pi(t)
and can be expressed as
(11)U(t)=∑i=02k‍liPi(t),
where
(12)li=-∑j=0k-1‍αi+1,jyn+j-∑j=1k‍αi+1,jyn+νj+hβi+1,kfn+k.By imposing the conditions (7) on (11), we obtain a system of 2k+1 equations, which can be expressed as
(13)QL=V,
where L=(l0,l1,…,l2k)T, V=(yn,yn+1,…,yn+k-1,yn+v1,…,yn+vk,fn+k)T, and Q is a matrix of dimension (2k+1)×(2k+1) defined as
(14)Q=(P0(tn)⋯P2k(tn)P0(tn+1)⋯P2k(tn+1)⋮⋮P0(tn+k-1)⋯P2k(tn+k-1)P0(tn+v1)⋯P2k(tn+v1)⋮⋮P0(tn+vk)⋯P2k(tn+vk)P0′(tn+k)⋯P2k′(tn+k)),
where Pj(t)∈Υm, j=0,…,2k, are basis functions. The elements of the vector L are obtained using Cramer’s Rule as
(15)lj=det⁡(Qj)det⁡(Q),j=0,1,…,2k,
where Qj is obtained by replacing the jth column of Q by V. The values of lj are substituted into (11) to obtain the continuous approximation
(16)U(t)=∑j=02k‍det⁡(Qj)det⁡(Q)Pj(t),
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 k=2 and for notational simplicity let {ν1,ν2}={r,s}, then on the interval from tn to tn+2=tn+2h, the approximation (6) becomes
(17)U(t)=-∑j=01‍αj(t)yn+j-∑j=12‍ανj(t)yn+νj+hβ2(t)fn+2,
with first derivative given by
(18)U′(t)=ddt(-∑j=01‍αj(t)yn+j-∑j=12‍ανj(t)yn+νj+hβ2(t)fn+2),
where αj(t), β2(t), and ανj(t) are continuous coefficients uniquely determined as
(19)α0=σ1(-12z+17rz-5r2z+17sz-24rszmm+7r2sz-5s2z+7rs2z-2r2s2z-16z2mm+11rz2-2r2z2+11sz2+rsz2-2r2sz2mm-2s2z2-2rs2z2+r2s2z2+4z3-9rz3mm+3r2z3-9sz3+6rsz3-r2sz3+3s2z3mm-rs2z3+8z4-3rz4-3sz4+rsz4),α1=σ2(20-28r+8r2-28s+39rs-11r2smm+8s2-11rs2+3r2s2+48z-48rz+12r2zmm-48sz+40rsz-9r2sz+12s2z-9rs2zmm+2r2s2z+24z2-8rz2-8sz2-6rsz2mm+3r2sz2+3rs2z2-r2s2z2-16z3+16rz3mm-4r2z3+16sz3-8rsz3+r2sz3mm-4s2z3+rs2z3-12z4+4rz4+4sz4-rsz4),αr=σ3(12z-17sz+5s2z+16z2-11sz2mm+2s2z2-4z3+9sz3-3s2z3-8z4+3sz4),αs=σ4(-12z+17rz-5r2z-16z2+11rz2mm-2r2z2+4z3-9rz3+3r2z3+8z4-3rz4),β2=ρ1-1(z-rz-sz+rsz+3z2-2rz2mmn-2sz2+rsz2+3z3-rz3-sz3+z4),
and where ρ1=20-8r-8s+3rs, σ1=(rsρ1)-1, σ2=[(r-1)(s-1)ρ1]-1, σ3=[(r-1)r(r-s)ρ1]-1, and σ4=[(r-s)(s-1)sρ1]-1. We also introduce the scale variable z for convenience only by letting z=(t-tn+1)/h.Suppose that the off-grid points are specified by the zeros of the second degree Chebyshev polynomial of the first kind given byr=(1/2)(2-2) and s=(1/2)(2+2); then (18) yields three additional methods when evaluated at t={tn+h,tn+hr,tn+hs} which are combined with the main method obtained by evaluating (17) at t=tn+2h to get a Chebyshev based method of the form (3) given as
(20)11hfn+1+hfn+2=-(4+92)yn+r-(4-92)yn+s+6yn+2yn+1,2(-2+2)hfn+2+22hfn+r=-2(12+52)yn+(8-262)yn+1-(-24+52)yn+r+(-8+112)yn+s,2(2+2)hfn+2-22hfn+s=-2(-12+52)yn+(8+112)yn+r-(24+52)yn+s+(8+262)yn+1,2hfn+2=-4(2+2)yn+s+yn+4yn+1+11yn+2+4(-2+2)yn+r,
whose stability properties are discussed in the next section.


## 2.1.k-Step HBDF

We define ak-step HBDF (main method) for the numerical solution of (1) on the interval [tn,tn+k] as
(2)∑j=0k‍αjyn+j+∑j=1k‍ανjyn+νj=hβkfn+k,
where k∈Z+ is the step number, tn=t0+nh, h=(tN-t0)/N, n=0,…,N, αj, ανj, βk are constant coefficients, and νj∈(0,k) are off-grid points which cannot be integers. We note that yn+j≈y(tn+j), yn+νj≈y(tn+νj), and fn+k=f(tn+k,yn+k). 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 k-step HBDF for solving (1).


## 2.2. Block HBDF

We define the blockk-step HBDF as
(3)∑j=0k-1‍α¯jyn+j+∑j=1k‍α¯νjyn+νj=hβ¯kfn+k+hβ¯ifn+i,i=1,…,k-1,∑j=0k-1‍α^jyn+j+∑j=1k‍α^νjyn+νj=hβ^kfn+k+hβ^νifn+νi,i=1,…,k,∑j=0k‍αjyn+j+∑j=1k‍ανjyn+νj=hβkfn+k,
where the first two members of (3) are additional methods provided by the continuous scheme constructed in the next subsection; αj, ανj, βνi, βk, α¯j, α¯νj, β¯νi, β¯k, α^j, α^νj, β^νi, β^k are constant coefficients, which are determined via the procedure given in the next subsection. The numerical solution yn+j is an approximation to the analytical solution y(tn+j), yn+νj≈y(tn+νj), fn+νj=f(tn+νj,yn+νj), and fn+j=f(tn+j,yn+j). In order to facilitate the analysis of (3) we express it in block form as follows:
(4)A1Wμ=A0Wμ-1+h[B1Fμ+B0Fμ-1],
where
(5)Wμ=(yn+ν1,…,yn+νk,yn+1,…,yn+k)T,Wμ-1=(yn-ν1,…,yn-νk,yn-k+1,…,yn)T,Fμ=(fn+ν1,…,fn+νk,fn+1,…,fn+k)T,Fμ-1=(fn-ν1,…,fn-νk,fn-k+1,…,fn)T,
where μ=1,…,Ϝ and n=0, k,…,N-k and Ai, Bi, i=0,1, are 2k by 2k 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ΥJ be the linear space of polynomials of degree ≤J and let U(t)∈ΥJ be a given polynomial. On the interval from tn to tn+k=tn+kh, the exact solution y(t) for (1) is approximated by
(6)U(t)=-∑j=0k-1‍αj(t)yn+j-∑j=1k‍ανj(t)yn+νj+hβk(t)fn+k,
where αvj(t) and βk(t) are continuous coefficients.Let the following conditions be satisfied:(7)U(tn+vj)=yn+vj,j=1,2,…,k,U(tn+j)=yn+j,j=0,1,2,…,k-1,U′(tk)=fn+k.
Further, suppose that (6) is defined by the polynomial basis functions
(8)αj(t)=∑i=02k‍αi+1,jPi(t),j={0,1,…,k-1},αvj(t)=∑i=02k‍αi+1,jPi(t),j={1,…,k},hβk(t)=∑i=02k‍hβi+1,kPi(t),
where the constants αi+1,j, αi+1,vj, and hβi+1,k are uniquely determined.Substituting (8) into (6) we get
(9)U(t)=-∑j=0k-1‍∑i=02k‍αi+1,jPi(t)yn+j-∑j=1k‍∑i=02k‍αi+1,jPi(t)yn+νj+h∑i=02k‍hβi+1,kPi(t)fn+k,
which simplifies to
(10)U(t)=∑i=02k‍{-∑j=0k-1‍αi+1,jyn+jmm-∑j=1k‍αi+1,jyn+νj+hβi+1,kfn+k}Pi(t)
and can be expressed as
(11)U(t)=∑i=02k‍liPi(t),
where
(12)li=-∑j=0k-1‍αi+1,jyn+j-∑j=1k‍αi+1,jyn+νj+hβi+1,kfn+k.By imposing the conditions (7) on (11), we obtain a system of 2k+1 equations, which can be expressed as
(13)QL=V,
where L=(l0,l1,…,l2k)T, V=(yn,yn+1,…,yn+k-1,yn+v1,…,yn+vk,fn+k)T, and Q is a matrix of dimension (2k+1)×(2k+1) defined as
(14)Q=(P0(tn)⋯P2k(tn)P0(tn+1)⋯P2k(tn+1)⋮⋮P0(tn+k-1)⋯P2k(tn+k-1)P0(tn+v1)⋯P2k(tn+v1)⋮⋮P0(tn+vk)⋯P2k(tn+vk)P0′(tn+k)⋯P2k′(tn+k)),
where Pj(t)∈Υm, j=0,…,2k, are basis functions. The elements of the vector L are obtained using Cramer’s Rule as
(15)lj=det⁡(Qj)det⁡(Q),j=0,1,…,2k,
where Qj is obtained by replacing the jth column of Q by V. The values of lj are substituted into (11) to obtain the continuous approximation
(16)U(t)=∑j=02k‍det⁡(Qj)det⁡(Q)Pj(t),
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 k=2 and for notational simplicity let {ν1,ν2}={r,s}, then on the interval from tn to tn+2=tn+2h, the approximation (6) becomes
(17)U(t)=-∑j=01‍αj(t)yn+j-∑j=12‍ανj(t)yn+νj+hβ2(t)fn+2,
with first derivative given by
(18)U′(t)=ddt(-∑j=01‍αj(t)yn+j-∑j=12‍ανj(t)yn+νj+hβ2(t)fn+2),
where αj(t), β2(t), and ανj(t) are continuous coefficients uniquely determined as
(19)α0=σ1(-12z+17rz-5r2z+17sz-24rszmm+7r2sz-5s2z+7rs2z-2r2s2z-16z2mm+11rz2-2r2z2+11sz2+rsz2-2r2sz2mm-2s2z2-2rs2z2+r2s2z2+4z3-9rz3mm+3r2z3-9sz3+6rsz3-r2sz3+3s2z3mm-rs2z3+8z4-3rz4-3sz4+rsz4),α1=σ2(20-28r+8r2-28s+39rs-11r2smm+8s2-11rs2+3r2s2+48z-48rz+12r2zmm-48sz+40rsz-9r2sz+12s2z-9rs2zmm+2r2s2z+24z2-8rz2-8sz2-6rsz2mm+3r2sz2+3rs2z2-r2s2z2-16z3+16rz3mm-4r2z3+16sz3-8rsz3+r2sz3mm-4s2z3+rs2z3-12z4+4rz4+4sz4-rsz4),αr=σ3(12z-17sz+5s2z+16z2-11sz2mm+2s2z2-4z3+9sz3-3s2z3-8z4+3sz4),αs=σ4(-12z+17rz-5r2z-16z2+11rz2mm-2r2z2+4z3-9rz3+3r2z3+8z4-3rz4),β2=ρ1-1(z-rz-sz+rsz+3z2-2rz2mmn-2sz2+rsz2+3z3-rz3-sz3+z4),
and where ρ1=20-8r-8s+3rs, σ1=(rsρ1)-1, σ2=[(r-1)(s-1)ρ1]-1, σ3=[(r-1)r(r-s)ρ1]-1, and σ4=[(r-s)(s-1)sρ1]-1. We also introduce the scale variable z for convenience only by letting z=(t-tn+1)/h.Suppose that the off-grid points are specified by the zeros of the second degree Chebyshev polynomial of the first kind given byr=(1/2)(2-2) and s=(1/2)(2+2); then (18) yields three additional methods when evaluated at t={tn+h,tn+hr,tn+hs} which are combined with the main method obtained by evaluating (17) at t=tn+2h to get a Chebyshev based method of the form (3) given as
(20)11hfn+1+hfn+2=-(4+92)yn+r-(4-92)yn+s+6yn+2yn+1,2(-2+2)hfn+2+22hfn+r=-2(12+52)yn+(8-262)yn+1-(-24+52)yn+r+(-8+112)yn+s,2(2+2)hfn+2-22hfn+s=-2(-12+52)yn+(8+112)yn+r-(24+52)yn+s+(8+262)yn+1,2hfn+2=-4(2+2)yn+s+yn+4yn+1+11yn+2+4(-2+2)yn+r,
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
(21)LTE=Zμ-(A1-1A0Zμ-1+hA1-1[B1F¯μ+B0F¯μ-1]),
where
(22)klZμ=(y(tn+ν1),…,y(tn+νk),y(tn+1),…,y(tn+k))T,klF¯μ=(f(tn+ν1,y(tn+ν1)),…,f(tn+νk,y(tn+νk)),f(tn+1,y(tn+1)),…,f(tn+k,y(tn+k)))T,Zμ-1=(y(tn-ν1),…,y(tn-νk),y(tn-k+1),…,y(tn))T,F¯μ-1=(f(tn-ν1,y(tn-ν1)),…,f(tn-νk,y(tn-νk)),f(tn-k+1,y(tn-k+1)),…,f(tn,y(tn)))T.
We recall here that the block method (4) is characterized by the matrices Ai, Bi, i=0,1, whose entries are given by the coefficients of (3).Definition 2.
The block method (4) has order of consistency p≥1.Thus, we expand the terms in (21) as a Taylor series about the point tn and also assume that the constant coefficients Cq=0 for q=0,1,…,p, leads to the local truncation error satisfying the maximum norm ∥LTE∥=Cp+1hp+1+O(hp+2).Theorem 3.
The block method (4) is zero-stable.Proof.
Zero-stability is concerned with the stability of a difference system in the limit ash tends to zero. Thus, as h→0 and noting that the first 2k-1 columns of A0 are all zeros with nonzero elements occurring only in the last column, the normalized first characteristic polynomial ρ(χ) of the method (4) is given by
(23)ρ(χ)=det⁡(χI-A1-1A0)=Γ1χ2K-1(χ+1).Γ1 is a constant and I is the identity matrix. Following Fatunla [15], the block method (4) is zero-stable, since from (23), ρ(χ)=0 satisfy |χj|≤1, j=1,…,2k, and for those roots with |χj|=1, 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 ρ(χ)=-(48/11)χ3(χ+1). Thus, the root condition is satisfied, since for ρ(χ)=0, χj={0,0,0,-1}, j=1,…,4, and for those roots with |χj|=1, 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 y′=λy gives the following discrete problem:
(24)Wμ=M(q)Wμ-1,q=λh,
where M(q)=(A1-qB1)-1(A0+qB0) is the stability matrix.Definition 5.
The block method (4) is said to be (i) A-stable if for all q∈C-, M(q) has a dominant eigenvalue qmax⁡ such that |qmax⁡|≤1; moreover, since qmax⁡ is a rational function, the real part of the zeros of qmax⁡ must be negative, while the real part of the poles of qmax⁡ must be positive; (ii) A0-stable if for all q∈R⊂C-, M(q) has a dominant eigenvalue qmax⁡ such that |qmax⁡|≤1; (iii) L0-stable if it is A0-stable and qmax⁡→0 as q→-∞; and (iv) L-stable if it is A-stable and qmax⁡→0 as q→-∞.Definition 6.
Let the stability region of (4) be Ω; then the method (4) is said to be A(α)-stable, α∈(0,π/2), if Ω⊃Ωα={q:|arg(-q)|<α,q≠0}.Theorem 7.
The block HBDF (4) is L0-stable and A(α)-stable.Proof.
The eigenvalues of the stability matrixM(q) in (24) are rational functions of q. It can be easily shown that the dominant eigenvalue qmax⁡ has the form qmax⁡=φ(q)/ψ(q), where the degree of the polynomial ψ(q) is greater than that of φ(q) and |φ(q)|≤|ψ(q)|. Applying Definitions 5 and 6 with α=90∘ completes the proof.Theorem 8.
The Chebyshev based method (20) is L0-stable and A(α)-stable.Proof.
The dominant eigenvalueqmax⁡ for the Chebyshev based method is given by
(25)qmax⁡=-48+36q+10q2+q348-60q+34q2-11q3+2q4.
The proof follows from Theorem 7.Moreover, the stability region for the Chebyshev based method is given in Figure1 showing the zeros and poles of the dominant eigenvalue qmax⁡.Figure 1
The region of absolute stability of the Chebyshev based HBDF of order 4 is to the left of the dividing line and is symmetric about the real axis; the square and plus symbols to the left and right of the imaginary axis represent the zeros and poles ofqmax⁡, respectively.


## 4. Numerical Application

### 4.1. Block Algorithm

Given the interval[t0,tN] where tn=t0+nh for n=0,1,…,N, the method (3) is implemented in block form (4) as follows.Step 1.
Letk=2; the method (3) reduces to the Chebyshev based method (20), and for μ=1, and n=0, the values of W1=(ur,u1,us,u2)T (T is the transpose) are simultaneously obtained over the subinterval [t0,t2] since W0 is known from the system of ODEs.Step 2.
Forμ=2, n=2, the values of W2=(ur+2,u3,us+2,u4)T are simultaneously obtained over the subinterval [t2,t4] since the values of W1 are known from the previous block.Step 3.
Forμ≥3, the process is continued to generate the vectors for the approximate solutions of the ODE on subintervals [t4,t6],…,[tN-2,tN].


### 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 2×2 and a 3×3 nonlinear stiff system.Example 1.
Consider the following linear system on the range0≤t≤1: y1′=-21y1+19y2-20y3, y2′=19y1-21y2+20y3, and y3′=40y1-40y2-40y3 with initial conditions y1(0)=1, y2(0)=0, and y3(0)=-1, respectively. The exact solution of the system is given by y1(t)=(1/2)(e-2t+e-40t(cos⁡⁡(40t)+sin⁡(40t))), y2(t)=(1/2)(e-2t-e-40t(cos⁡⁡(40t)+sin⁡(40t))), and y3(t)=(1/2)(2e-40t(sin⁡(40t)-cos⁡⁡(40t))).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 (p=4) 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 24=16. 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 (|y-y(t)|/(1+|y(t)|)).Table 1
A comparison of the Errs of methods for Example1.
Steps
HBDF
Rate
BVM
Rate
Ad-IVM
Rate
BDF
Rate
2
0
1.5
×
1
0
-
2
1.2
×
1
0
-
1
4.0
×
1
0
-
1
8.7
×
1
0
-
2
40
2.3
×
1
0
-
3
2.7
1.8
×
1
0
-
2
2.8
2.4
×
1
0
-
2
4.0
6.8
×
1
0
-
2
0.4
80
9.9
×
1
0
-
5
4.2
1.6
×
1
0
-
3
3.5
1.6
×
1
0
-
3
3.9
1.1
×
1
0
-
2
2.6
160
5.9
×
1
0
-
6
4.1
1.4
×
1
0
-
4
3.5
1.2
×
1
0
-
4
3.7
8.0
×
1
0
-
4
3.8
320
3.3
×
1
0
-
7
4.2
1.1
×
1
0
-
5
3.7
1.1
×
1
0
-
5
3.5
6.8
×
1
0
-
5
3.6
640
1.9
×
1
0
-
8
4.1
7.4
×
1
0
-
7
3.9
7.5
×
1
0
-
7
3.8
5.3
×
1
0
-
6
3.7Example 2.
We consider the following IVP which was solved by Vaquero and Vigo-Aguiar, [18]: y1′=-1002y1+1000y22 and y2′=y1-y2(1+y2) with initial conditions y1(0)=1 and y2(0)=1, respectively. The exact solution is y1(t)=e-2t; y2(t)=e-t.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].Table 2
A comparison of errors of methods att=5 for Example 2.
h
HBDF
Gauss-2s
EF-Gauss-2s
0.1
5.4907
×
1
0
-
11
5.1207
×
1
0
-
6
8.3555
×
1
0
-
7
0.01
6.0665
×
1
0
-
15
6.4579
×
1
0
-
12
4.2408
×
1
0
-
12Example 3.
We consider the given nonlinear system on the range0≤t≤48:
(26)y1′=-0.013y2-1000y1y2-2500y1y3,y1(0)=0,y2′=-0.013y2-1000y1y2,y2(0)=1,y3′=-2500y1y3,y3(0)=1.
True solution is
(27)yi(2)={-3.616933169289×10-6,mm0.9815029948230,1.018493388244-3.616933169289×10-6},yi(48)={-1.945338956808×10-6,-نن1.945338956808×10-6,kkk-3.616933169289×10-6,0.6110474831446,1.388950571516},mmmmmmmmmmmmmmmlmi=1,2,3
(see Jeltsch [19]).This example was chosen to show the performance of the HBDF on a3×3 nonlinear stiff system. In Table 3, the accuracy of the HBDF is measured by the maximum computing absolute errors at t={2,48}.Table 3
Absolute errors, (Err(yi)=|y(t)-yi|), i=1,2,3, for Example 3.
Method
h
t
=
2
t
=
48
Err(y1)
Err(y2)
Err(y3)
Err(y1)
Err(y2)
Err(y3)
HBDF
1
5.24
×
1
0
-
10
4.47
×
1
0
-
9
4.99
×
1
0
-
9
1.31
×
1
0
-
14
3.36
×
1
0
-
9
3.36
×
1
0
-
9
1/2
2.82
×
1
0
-
13
2.31
×
1
0
-
9
2.31
×
1
0
-
9
6.60
×
1
0
-
15
1.69
×
1
0
-
9
1.69
×
1
0
-
9
1/4
6.02
×
1
0
-
15
1.15
×
1
0
-
9
1.15
×
1
0
-
9
3.30
×
1
0
-
15
8.45
×
1
0
-
10
8.46
×
1
0
-
10
1/8
2.99
×
1
0
-
15
5.73
×
1
0
-
10
5.73
×
1
0
-
10
1.64
×
1
0
-
15
4.19
×
1
0
-
10
4.22
×
1
0
-
10
1/16
1.48
×
1
0
-
15
2.82
×
1
0
-
10
2.83
×
1
0
-
10
8.07
×
1
0
-
16
2.06
×
1
0
-
10
2.10
×
1
0
-
10


### 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
(28)∂u∂t=∂2u∂x2,(x,t)∈[0,1]×[0<t≤TN],
subject to the initial/boundary conditions u(x,0)=f(x), x∈[0,1], u(0,t)=u(1,t)=0, and t≥0.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 u(t)∈Rm be the exact solution of the resulting semidiscrete problem
(29)du(t)dt=Au(t),U(0)=ϕ,
where u(t)=[u1(t),u2(t),…,uM(t)]T, um(t)≈u(xm,t), ϕ=[ϕ(x1),ϕ(x2),…,ϕ(xM)]T, and
(30)A=1(Δx)2[-211-210⋱⋱⋱01-211-2].
The problem (29) is now a system of ordinary differential equations on which the block algorithm can be applied. We remark that if Δx is small, the system (29) is stiff.Example 4.
Consider the given PDE (see Cash [2])
(31)∂u∂t=κ∂2u∂x2,u(0,t)=u(1,t)=0,u(x,0)=sin⁡πx.
The exact solution u(x,t)=e-π2κtsin⁡πx.The output in Table4 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.Table 4
A comparison of errors of methods for Example4 at t=1.
Δ
x
Δ
t
κ
Crank-Nicolson
Cash (ab)
Cash (abc)
(HBDF)
0
.
1
0.1
1
3.0
×
1
0
-
5
1.5
×
1
0
-
5
4.5
×
1
0
-
6
1.43
×
1
0
-
6
0.1
0.1
5
2.0
×
1
0
-
4
3.0
×
1
0
-
8
2.0
×
1
0
-
10
3.66
×
1
0
-
15Example 5.
Consider the given stiff parabolic equation (see Cash [2])
(32)∂u∂t=κ∂2u∂x2,u(0,t)=u(1,t),u(x,0)=sin⁡πx+sin⁡ω¯πx,ω¯≫1.
The exact solution u(x,t)=e-π2κtsin⁡πx+e-ω¯2π2κtsin⁡ω¯πx.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 L0-stable are expected to perform poorly. The block HBDF is L-stable and hence performs well on this problem. Table 5 displays the results for κ=1 and a range of values for ω¯.Table 5
A comparison of errors of methods for Example5 at t=1 and ω-=1, Δx=0.1, Δt=0.1.
ω
-
(HBDF)
Crank-Nicolson
Cash (ab)
Cash (abc)
1
2.87
×
1
0
-
6
6.20
×
1
0
-
5
3.7
×
1
0
-
5
1.5
×
1
0
-
5
2
1.43
×
1
0
-
6
3.83
×
1
0
-
5
1.8
×
1
0
-
5
7.4
×
1
0
-
6
3
1.43
×
1
0
-
6
9.30
×
1
0
-
3
1.9
×
1
0
-
5
7.4
×
1
0
-
6
5
1.43
×
1
0
-
6
1.80
×
1
0
-
1
1.8
×
1
0
-
5
7.4
×
1
0
-
6
10
1.43
×
1
0
-
6
6.10
×
1
0
-
1
1.8
×
1
0
-
5
7.4
×
1
0
-
6


## 4.1. Block Algorithm

Given the interval[t0,tN] where tn=t0+nh for n=0,1,…,N, the method (3) is implemented in block form (4) as follows.Step 1.
Letk=2; the method (3) reduces to the Chebyshev based method (20), and for μ=1, and n=0, the values of W1=(ur,u1,us,u2)T (T is the transpose) are simultaneously obtained over the subinterval [t0,t2] since W0 is known from the system of ODEs.Step 2.
Forμ=2, n=2, the values of W2=(ur+2,u3,us+2,u4)T are simultaneously obtained over the subinterval [t2,t4] since the values of W1 are known from the previous block.Step 3.
Forμ≥3, the process is continued to generate the vectors for the approximate solutions of the ODE on subintervals [t4,t6],…,[tN-2,tN].


## 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 2×2 and a 3×3 nonlinear stiff system.Example 1.
Consider the following linear system on the range0≤t≤1: y1′=-21y1+19y2-20y3, y2′=19y1-21y2+20y3, and y3′=40y1-40y2-40y3 with initial conditions y1(0)=1, y2(0)=0, and y3(0)=-1, respectively. The exact solution of the system is given by y1(t)=(1/2)(e-2t+e-40t(cos⁡⁡(40t)+sin⁡(40t))), y2(t)=(1/2)(e-2t-e-40t(cos⁡⁡(40t)+sin⁡(40t))), and y3(t)=(1/2)(2e-40t(sin⁡(40t)-cos⁡⁡(40t))).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 (p=4) 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 24=16. 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 (|y-y(t)|/(1+|y(t)|)).Table 1
A comparison of the Errs of methods for Example1.
Steps
HBDF
Rate
BVM
Rate
Ad-IVM
Rate
BDF
Rate
2
0
1.5
×
1
0
-
2
1.2
×
1
0
-
1
4.0
×
1
0
-
1
8.7
×
1
0
-
2
40
2.3
×
1
0
-
3
2.7
1.8
×
1
0
-
2
2.8
2.4
×
1
0
-
2
4.0
6.8
×
1
0
-
2
0.4
80
9.9
×
1
0
-
5
4.2
1.6
×
1
0
-
3
3.5
1.6
×
1
0
-
3
3.9
1.1
×
1
0
-
2
2.6
160
5.9
×
1
0
-
6
4.1
1.4
×
1
0
-
4
3.5
1.2
×
1
0
-
4
3.7
8.0
×
1
0
-
4
3.8
320
3.3
×
1
0
-
7
4.2
1.1
×
1
0
-
5
3.7
1.1
×
1
0
-
5
3.5
6.8
×
1
0
-
5
3.6
640
1.9
×
1
0
-
8
4.1
7.4
×
1
0
-
7
3.9
7.5
×
1
0
-
7
3.8
5.3
×
1
0
-
6
3.7Example 2.
We consider the following IVP which was solved by Vaquero and Vigo-Aguiar, [18]: y1′=-1002y1+1000y22 and y2′=y1-y2(1+y2) with initial conditions y1(0)=1 and y2(0)=1, respectively. The exact solution is y1(t)=e-2t; y2(t)=e-t.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].Table 2
A comparison of errors of methods att=5 for Example 2.
h
HBDF
Gauss-2s
EF-Gauss-2s
0.1
5.4907
×
1
0
-
11
5.1207
×
1
0
-
6
8.3555
×
1
0
-
7
0.01
6.0665
×
1
0
-
15
6.4579
×
1
0
-
12
4.2408
×
1
0
-
12Example 3.
We consider the given nonlinear system on the range0≤t≤48:
(26)y1′=-0.013y2-1000y1y2-2500y1y3,y1(0)=0,y2′=-0.013y2-1000y1y2,y2(0)=1,y3′=-2500y1y3,y3(0)=1.
True solution is
(27)yi(2)={-3.616933169289×10-6,mm0.9815029948230,1.018493388244-3.616933169289×10-6},yi(48)={-1.945338956808×10-6,-نن1.945338956808×10-6,kkk-3.616933169289×10-6,0.6110474831446,1.388950571516},mmmmmmmmmmmmmmmlmi=1,2,3
(see Jeltsch [19]).This example was chosen to show the performance of the HBDF on a3×3 nonlinear stiff system. In Table 3, the accuracy of the HBDF is measured by the maximum computing absolute errors at t={2,48}.Table 3
Absolute errors, (Err(yi)=|y(t)-yi|), i=1,2,3, for Example 3.
Method
h
t
=
2
t
=
48
Err(y1)
Err(y2)
Err(y3)
Err(y1)
Err(y2)
Err(y3)
HBDF
1
5.24
×
1
0
-
10
4.47
×
1
0
-
9
4.99
×
1
0
-
9
1.31
×
1
0
-
14
3.36
×
1
0
-
9
3.36
×
1
0
-
9
1/2
2.82
×
1
0
-
13
2.31
×
1
0
-
9
2.31
×
1
0
-
9
6.60
×
1
0
-
15
1.69
×
1
0
-
9
1.69
×
1
0
-
9
1/4
6.02
×
1
0
-
15
1.15
×
1
0
-
9
1.15
×
1
0
-
9
3.30
×
1
0
-
15
8.45
×
1
0
-
10
8.46
×
1
0
-
10
1/8
2.99
×
1
0
-
15
5.73
×
1
0
-
10
5.73
×
1
0
-
10
1.64
×
1
0
-
15
4.19
×
1
0
-
10
4.22
×
1
0
-
10
1/16
1.48
×
1
0
-
15
2.82
×
1
0
-
10
2.83
×
1
0
-
10
8.07
×
1
0
-
16
2.06
×
1
0
-
10
2.10
×
1
0
-
10


## 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
(28)∂u∂t=∂2u∂x2,(x,t)∈[0,1]×[0<t≤TN],
subject to the initial/boundary conditions u(x,0)=f(x), x∈[0,1], u(0,t)=u(1,t)=0, and t≥0.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 u(t)∈Rm be the exact solution of the resulting semidiscrete problem
(29)du(t)dt=Au(t),U(0)=ϕ,
where u(t)=[u1(t),u2(t),…,uM(t)]T, um(t)≈u(xm,t), ϕ=[ϕ(x1),ϕ(x2),…,ϕ(xM)]T, and
(30)A=1(Δx)2[-211-210⋱⋱⋱01-211-2].
The problem (29) is now a system of ordinary differential equations on which the block algorithm can be applied. We remark that if Δx is small, the system (29) is stiff.Example 4.
Consider the given PDE (see Cash [2])
(31)∂u∂t=κ∂2u∂x2,u(0,t)=u(1,t)=0,u(x,0)=sin⁡πx.
The exact solution u(x,t)=e-π2κtsin⁡πx.The output in Table4 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.Table 4
A comparison of errors of methods for Example4 at t=1.
Δ
x
Δ
t
κ
Crank-Nicolson
Cash (ab)
Cash (abc)
(HBDF)
0
.
1
0.1
1
3.0
×
1
0
-
5
1.5
×
1
0
-
5
4.5
×
1
0
-
6
1.43
×
1
0
-
6
0.1
0.1
5
2.0
×
1
0
-
4
3.0
×
1
0
-
8
2.0
×
1
0
-
10
3.66
×
1
0
-
15Example 5.
Consider the given stiff parabolic equation (see Cash [2])
(32)∂u∂t=κ∂2u∂x2,u(0,t)=u(1,t),u(x,0)=sin⁡πx+sin⁡ω¯πx,ω¯≫1.
The exact solution u(x,t)=e-π2κtsin⁡πx+e-ω¯2π2κtsin⁡ω¯πx.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 L0-stable are expected to perform poorly. The block HBDF is L-stable and hence performs well on this problem. Table 5 displays the results for κ=1 and a range of values for ω¯.Table 5
A comparison of errors of methods for Example5 at t=1 and ω-=1, Δx=0.1, Δt=0.1.
ω
-
(HBDF)
Crank-Nicolson
Cash (ab)
Cash (abc)
1
2.87
×
1
0
-
6
6.20
×
1
0
-
5
3.7
×
1
0
-
5
1.5
×
1
0
-
5
2
1.43
×
1
0
-
6
3.83
×
1
0
-
5
1.8
×
1
0
-
5
7.4
×
1
0
-
6
3
1.43
×
1
0
-
6
9.30
×
1
0
-
3
1.9
×
1
0
-
5
7.4
×
1
0
-
6
5
1.43
×
1
0
-
6
1.80
×
1
0
-
1
1.8
×
1
0
-
5
7.4
×
1
0
-
6
10
1.43
×
1
0
-
6
6.10
×
1
0
-
1
1.8
×
1
0
-
5
7.4
×
1
0
-
6


## 5. Conclusion

In this paper, we have proposed a generalized high order block hybridk-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.



---

*Source: 162103-2014-10-20.xml*