Printer Friendly


Byline: M. A. Rehmana and M. S. A. Taj


This paper deals with numerical method for the approximate solution of one-dimensional heat equation (Equation) with integral boundary conditions. The integral conditions are approximated by using Simpson's 1/3 rule while the space derivatives are approximated by third-order finite difference approximations. Then method of lines, semidiscritization approach, is used to transform the model partial differential equation into a system of first-order linear ordinary differential equations whose solution satisfies a recurrence relation involving matrix exponential function. The method developed is L-acceptable, third-order accurate in space and time and do not require the use of complex arithmetic. A parallel algorithm is also developed and implemented on several problems from literature and found to be highly accurate when compared with the exact ones and alternative techniques.

Keywords: Heat equation, Boundary integral specifications, Third order numerical methods, Method of lines, Parallel algorithm


In this paper we have considered the one-dimension heat equation with non-local boundary conditions. Much attention has been paid in the literature for the development, analysis and implementation of accurate methods for the numerical solution of this problem. Consider heat equation


where (Equation) and (Equation) are known functions and are assumed to be sufficiently smooth to produce a smooth classical solution of(Equation). is given positive constant. A number of numerical procedures are suggested in the literature to address the problem: see, for instant [2, 3, 5, 6, 9, 12].

Inspiring from great accuracy achieved in [11] the authors aim to attempt this problem. In this paper the method of lines, semi discretization approach, will be used to transform the model partial differential equation (PDE) into a system of first-order, linear, ordinary differential equations (ODEs), the solution of which satisfies a recurrence relation involving matrix exponential terms. A third-order rational approximation will be used to approximate exponential functions which will lead to an -acceptable algorithm which may be parallelized through the partial fraction splitting technique.


Dividing the interval(Equation) into (Equation) subintervals each of width (Equation), so that (Equation) and the time variable (Equation) into time steps each of length (Equation) gives a rectangular mesh of points with co-ordinates (Equation) where (m = 0, 1, 2, ..., N+1 and n = 0, 1, 2, ...) covering the region (Equation) and its boundary (Equation) consisting of lines x = 0, x = 1 and t=0.

The space derivative in (1) may be approximated to the third-order accuracy at some general point (x, t) of the mesh by using five point central difference approximations.


It is worth noting that the equation (5) is valid only for (Equation) with (Equation) To attain the same accuracy at the points (Equation) for (Equation) and (Equation), special formulae developed by [11] are used, which approximate (Equation) not only to third-order but also with dominant error term (Equation) for (Equation) and (Equation). Such approximations are



for the mesh points (Equation) and (Equation) respectively. Applying (1) with (5), (6) and (7) to all the interior mesh points of the grid at time level (Equation) produces a system of N linear equations in N+2 unknowns . The integrals in (3) and (4) are approximated by using Simpson's (Equation) rule as used by [6, 7, 8]. Here

(Equation) (9)

where N is an odd integer. Solving (8) and (9) for U0 and UN+1 and substituting these values in the above system we have a system of N linear ordinary differential equations with N unknowns U1, ..., UN which can be written in vector matrix form as

(Equation) (10)

with initial distribution

(Equation) (11)

in which (Equation), where T denoting the transpose and matrix A of order NxN which is given by

also (Equation) and (Equation). The column vector v(t) contains the contribution from the functions (Equation) and is given as


The solution of the system (10) subject to (11) is given by

(Equation) (13)

which satisfies the recurrence relation

(Equation) (14)

Eigen-values of the matrix A are calculated using MATLAB for N=9, 19, 39, 79 and it is observed that they are distinct negative real ones or complex with negative real parts.

To approximate the matrix exponential in (14) following [11] we use a rational approximation given by

(Equation) (15)

(Equation) is a real scalar. The coefficients and (Equation) and (Equation) are real and interlinked as (Equation) and (Equation). So we have

(Equation) (16)


(Equation) (17)

The quadrature term in (14) is approximated as

(Equation) (18)

where (Equation) and W1, W2, W3 are matrices given by [11] and are mentioned here for convenience.

(Equation) (19)

(Equation) (20)

(Equation) (21)


(Equation) (22)


In this section the numerical method will be applied to a problem from the literature and results obtained will be compared with exact solution as well as with the results existing in the literature. Following [11] we have chosen here (Equation) and (Equation) as produced in [10] which gives (Equation). It is found that r1 =2.35913888475789, r2 = 2.34684670205017 and r3 =2.17953792391154 are the real zeros of (Equation) for this choice. For this choice of values L-stability is also guaranteed.

Example (1):- Consider the problem (1)-(4) with


which has theoretical solution (Equation)

The problem is solved using the scheme developed in this paper for h = l = 0.05, 0.025, 0.01, 0.005, 0.0025, 0.001 at x = 0.6 and t =1. The relative errors obtained by this scheme are given in Table 1 and the results are compared with different schemes, BTCS implicit scheme, Crandall method, FTCS scheme and Dufort-Frankel scheme given by [7]. From the table we can see that the results of the new scheme are far better than those of the schemes given in [7].

Example (2):- Now consider the problem (1)-(4) with


where (Equation) and (Equation)

which has theoretical solution

(Equation) [7].

For Example (2) results are given in Table 2 and Table 3. In Table 2 the results are computed for h = l = 0.05, 0.025, 0.01, 0.005, 0.0025, 0.001 at x = 0.6 and t = 0.1. The relative errors developed in the scheme are compared with different schemes, BTCS implicit scheme, Crandall method, FTCS scheme and Dufort-Frankel scheme given by [7]. From the table it is clear that the results are in good agreement as compared with the exact ones as well as better than other schemes. Moreover the new scheme is third-order accurate except for very small values of h and l where accumulating error is high.

The example (2) is also solved for h=l=0.01 for different values of t at x=0.25 and the results are presented in Table 3. Table 3 shows that the scheme developed in this paper gives superior results to those computed by using the Crank-Nicolson finite-difference method [12], the implicit finite-difference technique and the parallel techniques [6]. The parallel technique developed in [6] is second-order accurate while the parallel technique developed in this paper is third-order accurate.

Example (3):- Again consider the problem (1)-(4) with


In Example (3) results computed are given in Table 4 and Table 5. In Table 4 results are calculated for h = 0.01 = l at x = 1 and for different values of t. From the table it is clear that the analytical solution calculated by using the scheme developed in this paper is good agreement with the exact ones. Also the solution converges towards exact solution as t increases.

In Table 5 results are given for t = 1 with h = l = 0.1, 0.05, 0.025, 0.0125 and 0.00625 at x = 0.5 and x = 1. It is clear from the table that results for x = 1 are far better than those at x = 0.5 also the method is third-order accurate. CPU time taken for the new scheme developed in this paper is also given in the table which shows that the new scheme is very fast.

This problem is also solved by [2] with (Equation) and (Equation). It is noted that the increase in the value of causes more accuracy.


It is observed that the results obtained using new scheme are highly accurate as compared to those of other schemes and the method developed is third-order accurate in space and time as well as L-acceptable. This technique can be coded easily on serial or parallel computers.

It is worth mentioning that the method using real arithmetic and multiprocessor architecture especially in multi-dimensional problems will save remarkable CPU time rather than the complex arithmetic based methods.


[1] M. Akram, and M. A. Pasha, " A numerical method for the heat equation with a nonlocal boundary condition", Inter. J. Information and systems sciences, 1(2),162-171(2005).

[2] W. T. Ang. "A method of solution for the one-dimensional heat equation subject to nonlocal conditions".SEA Bull. Math. 26(2),197-203(2002).

[3] N Borovykh. "Stability in the numerical solution of the heat equation with nonlocal boundary conditions". Appl. Numer. Math.,4217-27(2002).

[4] W. A. Day, " Extension of a property of the heat equation to linear thermoelasticity and other theories", Quart. Appl. Math., 40, 319-330(1982).

[5] M. Deghan, "On the numerical solution of diffusion equation with a non-local boundary condition", Math. Prob. in Engg., 2, 81-92(2003).

[6] M. Deghan," Numerical solution of a parabolic equation with non-local boundary specifications". Appl. Math. Comp.,145, 185-194(2003).

[7] M. Deghan," Efficient techniques for the second-order parabolic equation subject to nonlocal specifications", Appl. Numer. Math., 52, 39-62(2005).

[8] A. B. Gumel,"On the numerical solution of the diffusion equation subject to the specification of mass". J. Austral. Math. Soc. Ser., B 40, 475-483(1999).

[9] Y. Liu, "Numerical solution of the heat equation with nonlocal boundary conditions". J. Comp. Appl. Math.,110, 115-127(1999).

[10] M. S. A. Taj, PhD. Thesis, Brunel University Uxbridge, Middlesex, England (1995).

[11] M. S. A. Taj and E. H. Twizell, "A family of third-order parallel splitting methods for prabolic partial differential equations", Intern. J. Comput. Math., 67, 411-433(1998).

[12] S. Wang, and Y. Lin." A numerical method for the diffusion equation with nonlocal boundary specifications". Inter. J. Engrg. Sci., 28(1990)543-546.

Department of Mathematics, University of Management and Technology, Lahore, Pakistan. a, Department of Mathematics, COMSATS Institute of information Technology, Lahore, Pakistan. b, a, b
COPYRIGHT 2012 Asianet-Pakistan
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2012 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Rehmana, M. A.; Taj, M. S. A.
Publication:Science International
Article Type:Report
Geographic Code:9PAKI
Date:Mar 31, 2012

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters