L-Stable Block Hybrid Second Derivative Algorithm for Parabolic Partial Differential Equations

An L-stable block method based on hybrid second derivative algorithm (BHSDA) is provided by a continuous second derivative method that is defined for all values of the independent variable and applied to parabolic partial differential equations (PDEs). The use of the BHSDA to solve PDEs is facilitated by the method of lines which involves making an approximation to the space derivatives, and hence reducing the problem to that of solving a time-dependent system of first order initial value ordinary differential equations. The stability properties of the method is examined and some numerical results presented.

We seek a solution in the strip [ ] [ ] 0,1 0 < t T × ≤ by first fixing the grid in the spatial variable x , then approximating this spatial derivative using the central difference method, and finally solving the resulting system of first order time dependent ODEs.Specifically, we discretize the space variable with mesh spacings We then define , and replace the partial derivatives ( ) , , , which reduces the PDE to the semidiscrete problem ( ) ( ) which can be written in the form , t = f u Au , and A is an M M × matrix arising from the central difference approximations to the derivatives of x .The problem ( 2) is now a system of first order ODEs which is solved by the BHSDA.
The paper is organized as follows.In Section 2, we derive a continuous approximation which is used to obtain the BHSDA.The BHSDA is also analyzed in Section 2. The computational aspects of the method is given in Section 3. Numerical examples are given in Section 4 to show the accuracy of the method.Finally, the conclusion of the paper is discussed in Section 5.

Development of the Method
We begin by considering a scalar form of (3) where we assume that the function f is Lipshitz continuous and the problem (4) possesses a unique solution.Furthermore, let n u be an approximation of the theoretical solution ( ) where j  are unknown coefficients.In order to uniquely determine the unknown coefficients j  , we impose that the interpolating function (4) coincides with the analytical solution at the end point n t and also satisfies the differential Equation (3) at the points , 0,1, 2 to obtain the following system of equations: , , , 0,1, 2.
We note that (6) leads to a system of five equations which is solved by Cramer's Rule to obtain j  .The continuous method is constructed by substituting the values of j  into Equation ( 5) which is simplified and expressed in the form where ( ) Remark 2.1 In order to conveniently analyze and implement the method (8), we will express it in block form as given in (9). where  , and the matrices ( ) are 2 by 2 matrices whose entries are given by the coefficients of (8).

Local Truncation Error
Define the local truncation error of (4) as where ( ) , , , , , , is a linear difference operator.Assuming that ( ) z t is sufficiently differentiable, we can expand the terms in (10) as a Taylor series about the point n t to obtain the expression for the local truncation error.
( ) ( ) , hence the method is of order four.

Stability
Proposition 2. 2 The BHSDA (9) applied to the test equations u u with the amplification matrix ( ) Remark 2. 3 The dominant eigenvalue of ) (q M specified by 2 max 2 3 48 18 2 48 30 8 q q q q q q + + = − + − is a rational function called the stability function which determines the stability of the method. Proof.We begin by applying (2) to the test equations u u λ ′ = and , we obtain a system of linear equations which is used to solve for Y µ with (12) as a consequence.

Definition 2.4
The block method ( 9) is said to be 1) A -stable if for all q − ∈ ,

( )
M q has a dominant eigenvalue max q such that max 1 q ≤ ; moreover, since max q is a rational function, the real part of the zeros of max q must be negative, while the real part of the poles of max q must be positive; 2) L -stable if it is A -stable and max 0 q → as q → −∞ .Corollary 2. 5 The method ( 9) is A -stable and L -stable.
Proof: The dominant eigenvalue max q for the method ( 9) is given by 48 18 2 48 30 8 q q q q q q + + = − + − and the proof follows from definition 2.4.Remark 2. 6 The stability region for the method ( 9) is given in Figure 1 showing the zeros and poles of the dominant eigenvalue max q .

Computational Aspects
The resulting system of ODEs ( 3) is then solved on the partition  , N is a positive integer and n the grid index.
Step 1: Use the block method (9) to solve (3 Step , Step 3:  , to generate the approxi- mations , We note that for linear problems, we solve (3) directly with our Mathematica code enhanced by the feature

Numerical Examples
Computations were carried out in Mathematica 9.0 and the errors were calculated as , where . We note that the method is particularly useful, but not limited to solving parabolic partial differential equations where the solution decays very rapidly and where the PDEs are stiff parabolic equations (see Cash [4]).
Example 4.1 As our first test example, we solve the given PDE (see Cash [4]) The exact solution ( ) In Table 1, it is noticed that the method with the BHSDA is the most accurate.Example 4.2 As our second test example, we solve the given stiff parabolic equation (see Cash [4]) 2 , 0, 1, = 0, , 0 sin π sin π , 1.
The exact solution ( ) Cash [4] notes that as ω increases, equations of the type given in example 4.2 exhibit characteristics similar to model stiff equations.Hence, the methods such as the Crank-Nicolson method which are not 0 L -stable are expected to perform poorly.The BHSDA is L -stable and perform excellently when applied to this problem.Therefore the BHSDA is competitive with the 0 L -stable methods of Cash [4].In Table 2, we display the results for 1 κ = and a range of values for ω .

Conclusion
We have proposed a BHSDA for solving parabolic PDEs via the method of lines.The method is shown to be L - stable and competitive with existing methods in the literature.
1) by the central difference approximation to obtain q hλ =

Figure 1 .
Figure1.The region of absolute stability of the BHSDA 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 of q max respectively.

Table 1 .
A comparison of errors of methods for Example 4.1 at t = 1.