ON THE NUMERICAL SOLUTION OF THE DIFFUSION EQUATION WITH A NONLOCAL BOUNDARY CONDITION

Parabolic partial differential equations with nonlocal boundary specifications feature in the mathematical modeling of many phenomena. In this paper, numerical schemes are developed for obtaining approximate solutions to the initial boundary value problem for one-dimensional diffusion equation with a nonlocal constraint in place of one of the standard boundary conditions. The method of lines (MOL) semidiscretization approach is used to transform the model partial differential equation into a system of first-order linear ordinary differential equations (ODEs). The partial derivative with respect to the space variable is approximated by a second-order finite-difference approximation. The solution of the resulting system of first-order ODEs satisfies a recurrence relation which involves a matrix exponential function. Numerical techniques are developed by approximating the exponential matrix function in this recurrence relation. We use a partial fraction expansion to compute the matrix exponential function via Pade approximations, which is particularly useful in parallel processing. The algorithm is tested on a model problem from the literature.


Introduction
Over the last few years, many physical phenomena were formulated into nonlocal mathematical models [1,2,3,4,5,6,7,8].However, such problems were not well studied in general.For instance, there were few discussions on the numerical solutions of the parabolic partial differential equations with nonlocal boundary specifications.
In this paper, we are concerned with one-dimensional parabolic equations with a nonlocal condition: the so-called energy specification.This is a linear constraint having the form b 0 u(x,t)dx = E(t), where b is a constant and E is a given function.Coupled with a one-dimensional parabolic equation, this condition is quite different from the classical boundary condition.
The nonlocal problems are very important in the transport of reactive and passive contaminates in aquifer, an area of active interdisciplinary research of mathematicians, engineers, and life scientists.We refer the reader to [9,10] for the derivation of mathematical models and for the precise hypothesis and analysis.
Mathematical formulation of this kind also arises naturally in various engineering models, such as nonlocal reactive transport in underground water flows in porous media [9,11], heat conduction, radioactive nuclear decay in fluid flows [25], non-Newtonian fluid flows, viscoelastic deformations of materials with memory (in particular polymers) [24], semiconductor modeling [1], and biotechnology.
The number of physical phenomena modeled by partial differential equations involving nonlocal integral terms is constantly increasing.The authors of [23] have given an example from metrology.This example is a model for the evolution of the temperature distribution of air near the ground during calm clear nights.
Parabolic problems with nonlocal boundary conditions also arise in the quasistatic theory of thermoelasticity [12,13].An interesting collection of nonlocal parabolic problems in one-space dimension is discussed in [19].
One very important characteristic of all these models is that they all express a conservation of a certain quantity (mass, momentum, heat, etc.) in any moment for any subdomain.This in many applications is the most desirable feature of the approximation method when it comes to the numerical solution of the corresponding initial boundary value problem.
Only recently much attention has been paid in the literature to the development, analysis, and implementation of accurate methods for the numerical solution of timedependent partial differential equations with nonlocal boundary conditions [14,15,16].
The familiar finite-difference, finite-element methods have been used by many authors in the last ten years for this model [4] and some other similar problems [2,3,5,17,18,28].
The purpose of this paper is to use the Pade approximant for solving the one-dimensional time-dependent diffusion equation subjected to the given initial condition the boundary condition and the nonlocal boundary condition where f , g, b, s, and m are known, while the function u is to be determined.The existence, uniqueness, and continuous dependence on the data of the solution of this nonlocal problem has been studied in [6,20].

Mehdi Dehghan 83
This paper is organized in the following way.Numerical schemes for the solution of (1.1), (1.2), (1.3), and (1.4) are described in Section 2. Development of numerical methods based on the Pade approximants is presented in Section 3. A discussion on the parallel algorithms is given in Section 4. To approximate the matrix exponential function, a rational approximation consisting of three parameters is described in this section.The numerical results for two test problems produced by the methods developed are given in Section 5. Section 6 concludes this paper with a brief summary.

Treatment of the nonlocal boundary specification
We reformulate the nonlocal boundary condition (1.4) to the separated Neumann type condition first by differentiating with respect to t, and then using (1.1).Thus (2.1) serves as the boundary condition at zero.Equations (1.1), (1.2), (1.3), and (2.1) can now be spatially differenced in the standard way, provided that b is a mesh point.We divide the domain [0,1] × [0,T] into an M × N mesh with spatial step size h = 1/M in x direction and the time step size l = T/N, respectively.Grid points (x i ,t n ) are given by in which M is an integer.We use u n i to denote the approximation of u(ih,nl).The solution vector U n will be denoted in the following form: In this paper, the method of lines semidiscretization approach is used to transform the model partial differential equation into a system of first-order linear ordinary differential equations, the solution of which satisfies a certain recurrence relation involving matrix exponential terms.The development of numerical methods will be based on Pade approximations to such exponentials.A technique which does not require the use of complex arithmetic is used to approximate the solution of the one-dimensional heat equation at interior grid points.
The space derivative in (1.1) and boundary terms will be replaced by the finite-difference approximation It is worth noting that (2.4) is useful for i = 1,2,...,M − 2. To attain the same accuracy at the points (x i ,t n ) with i = M − 1 and i = M, the approximations (2.5) Applying (1.1) to all the (M − 1) interior mesh points of the interval [0,1], at time level t n = nl, with the space derivative replaced by (2.4), produces a system of (M − 1) first-order, linear ordinary differential equations of the form with the initial condition in which the matrix A is in the following form: (2.8) Note that the inhomogeneous term ψ(t) contains the contribution from the boundary functions m (t) and g(t).
Solving the system of ODEs (2.6) subject to the initial condition (2.7) gives which satisfies [22] the recurrence relation in which l is a constant time step in the discretization of the time variable t ≥ 0 at the points t n = nl (n = 0,1,2,...).

Solution at the first time step and the Pade approximant
To approximate the matrix exponential in (2.10), we consider rational approximation to e θ of the form where the coefficients a i and b i are real and and where Note that the denominator of R(θ) has distinct real zeros provided that a 2 2 − 3a 1 a 3 > 0. We have chosen a 1 = 2, a 2 = 14/15, and a 3 = 1/10.It can be shown that using these values, L-stability is granted [27].
The quadrature term in (2.10) will be approximated by where s 1 = s 2 = s 3 and W 1 , W 2 , and W 3 are matrices.
Assuming that r i (i = 1,2,3) are distinct real zeros of q(θ), then and the approximation to the matrix exponential function may be written in partial fraction form as in which p i , the partial fraction coefficients of R(θ), are defined by Note that the implementation of the method using a parallel architecture with three processors is based on the partial fraction decomposition (see [21]) of exp(lA)u(t), W 1 ψ(t), W 2 ψ(t + l/2), and W 3 ψ(t + l) in (3.12).
(3.18) So, U(t + l) = y 1 + y 2 + y 3 in which y 1 , y 2 , and y 3 are, respectively, the solutions of the systems

The parallel algorithms
Note that (3.19) have great importance in the parallel environment since they can be used to solve the corresponding linear algebraic systems on processors operating concurrently.U(t + l) in (3.12), the solution vector at time t = (n + 1)l, may now be obtained via the parallel algorithms using three different processors in the following form.
Step 5: Use Step 6: Solve Then Mehdi Dehghan 89  In implementing these algorithms, Processor I generates once-only decomposition of I − (l/r 1 )A, while Processor II generates once-only decomposition of I − (l/r 2 )A and Processor III generates once-only decomposition of I − (l/r 3 )A.

Numerical test
The numerical method described in the previous sections was applied to two problems from the literature.
which is easily seen to have exact solution (5. 2) The results of u(0.25,1.0)with h=0.01,0.005,0.0025,0.001,and l = 0.01,0.005,0.0025,0.001, using the scheme discussed in Section 4, are shown in Table 5.1 and are compared with the results obtained using the implicit finite-difference scheme of [4].
In Table 5.1, we present the relative error, |Uapprox − U exact |/|U exact |, at x = 0.25 for our formula and the method of [4].
The results obtained using the new scheme developed in this paper are slightly more accurate than those from the scheme of [4].Note that the new scheme will require less CPU time.It is clear that as far as efficiency is concerned, the scheme introduced in this which is easily seen to have the exact solution u(x,t) = exp − π 2 t sin(πx). (5.4) The results for our second example are given in Table 5.2.Calculations were performed for h = l = 0.01,0.005,0.0025,0.001.
Overall, the results showed that the technique developed in this paper gave superior numerical results to those computed using the method of [4].It is worth mentioning that the use of only real arithmetic especially in multispace dimensions can yield large saving in CPU time used.

Conclusion
In this paper, an algorithm was applied to the one-dimensional diffusion equation with a nonlocal condition replacing one standard boundary condition.The second-order spatial derivative was discretized to result in an approximating system of ODEs.The exact solution of this system of first-order ODEs satisfies a recurrence relation involving the matrix exponential function.This function is approximated by a rational function possessing real and distinct poles, which consequently readily admits a partial fraction expansion, thereby allowing the distribution of the work in solving the corresponding linear algebraic Mehdi Dehghan 91 systems on concurrent processors.The method developed does not require the use of complex arithmetic and needs only real arithmetic in its implementation.This technique worked very well for one-dimensional diffusion with an integral condition.For the model problem considered, the parallel algorithm (4.7) was found to be about three times faster than the standard implicit finite-difference scheme of [4].A comparison with the standard implicit scheme for the model problem clearly demonstrates that the new technique is computationally superior.The numerical test obtained by using the method described in this paper gives acceptable results and suggests convergence to the exact solution when h goes to zero.
2.0 × 10 −6 paper is a better candidate for the model problem.This technique can be coded very efficiently on super/parallel computers.