International Journal of Computational Mathematics
Volume 2014 (2014), Article ID 923693, 9 pages
http://dx.doi.org/10.1155/2014/923693

Research Article

A Numerical Method for 1-D Parabolic Equation with Nonlocal Boundary Conditions

M. Tadi¹ and Miloje Radenkovic²

¹Department of Mechanical Engineering, University of Colorado at Denver, Campus Box 112, P.O. Box 173364, Denver, CO 80217-3364, USA
²Department of Electrical Engineering, University of Colorado at Denver, Campus Box 112, P.O. Box 173364, Denver, CO 80217-3364, USA

Received 11 July 2014; Revised 2 November 2014; Accepted 3 November 2014; Published 20 November 2014

Academic Editor: Chengpeng Bi

Copyright © 2014 M. Tadi and Miloje Radenkovic. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper is concerned with a local method for the solution of one-dimensional parabolic equation with nonlocal boundary conditions. The method uses a coordinate transformation. After the coordinate transformation, it is then possible to obtain exact solutions for the resulting equations in terms of the local variables. These exact solutions are in terms of constants of integration that are unknown. By imposing the given boundary conditions and smoothness requirements for the solution, it is possible to furnish a set of linearly independent conditions that can be used to solve for the constants of integration. A number of examples are used to study the applicability of the method. In particular, three nonlinear problems are used to show the novelty of the method.

1. Introduction

In this paper we consider a numerical method for 1-D heat conduction problems with nonlocal boundary conditions. Such problems appear very naturally in a number of physical systems including heat conduction [1], elastic deformation of thermoelastic rods [2, 3], chemical reactions [4, 5], population dynamics [6], and petroleum exploration [7]. A recent review of problems with nonlocal boundary conditions can also be found in [8] and references therein. Traditional finite difference methods have been shown to have problems in accuracy [9]; therefore, a number of investigators have developed various numerical methods for the above problem. Recent results include methods based on -based finite difference [10], reproducing kernel space [11], and Adomian expansion [12].

The purpose of this paper is to apply a method based on local coordinates. We have developed this method for the solution of a Stefan problem [13]. The algorithm transforms the working equations onto a local coordinate system. It is then possible to write down an exact solution which is valid and can be used within a small local region. The procedure leads to an implicit scheme that is first order accurate in time. However, the accuracy in time can be easily improved. The novelty of the method is in the fact that it can obtain exact solutions in space based on local coordinates. In addition, the algorithm can be applied to a large class of nonlinear problems. The formulation provides a natural way to obtain a numerical solution to the problem without any iterations which is often the case using the existing methods. The present method can provide an effective way to numerically investigate the behavior of linear and nonlinear singular parabolic equations with nonlocal source terms that lead to blowup [14, 15].

Section 2 introduces the method in detail for a 1-D parabolic problem. Section 3 studies the stability of the method and the invertibility of the coefficient matrix. Section 4 uses a number of examples to investigate the applicability of the method and compares the results to exact solutions. In particular, three nonlinear problems are studied. The last two nonlinear problems are used to show the novelty of the present method.

2. A numerical Method Based on Local Coordinates

Consider a one-dimensional heat equation with nonlocal boundary conditions given by

where the term is the heat generation, is the initial condition, and the boundary conditions are of nonlocal type. Assume that the domain is divided into equal intervals. Consider a finite difference mesh given in Figure 1.

Figure 1: Finite difference mesh showing the local coordinate.

For a typical node , consider an implicit expansion of the equation at time interval given by

where the value of is known. Now, consider a transformation to the local coordinate, , according to

In terms of the local coordinate, (2) simplifies to

The variable denotes the heat generation function evaluated at time () and location . (In general, it is a function of space ; however, for a first-order accurate method, the value at can be used. Using Taylor series expansion, it is possible to obtain a more accurate description in terms of ; that is, , for which exact solution can be obtained.) Treating as a local constant in function allows us to write down the exact solution of around every node given by

where now and are constants of integration. For intervals, there are nodes. Therefore, there are equations similar to (5), which describe the temperature profile; that is, . It follows that there are constants, that is, () that are unknown and should be computed at each time interval. It is possible to introduce linearly independent conditions. These conditions are not unique. However, a linearly independent set of conditions is given by the following.

(i)The temperature distribution at should be continuous. We can impose this condition at mid-points between nodes. It leads to

(ii)The flux at should be continuous. In the present case, it follows that the temperature at should have a continuous first derivative. We can also impose this condition at mid-points between nodes. It leads to

The above conditions furnish linearly independent equations. The additional two conditions needed are given by the boundary conditions. The solution must satisfy the boundary conditions at and . These conditions are nonlocal and can be imposed according to

The integral terms on the right hand side of the above equations can be computed in terms of the solutions given in (5) according to the following convention. For the first node, the domain of integration can be chosen to be . For the nodes , the domains can be chosen to be . For the last node, the domain can be chosen as . In terms of the local coordinates , for the interior nodes, the individual domains of integrations are all . Similarly, the domain of integration is for the first node and for the last node. Using these domains, the integral condition at is given by

and for , it is given by

The functions and are known, and the above integrals can be evaluated in closed forms. After evaluating the integrals, (10) leads to

Note that is known for all nodes, and the integral on the right hand side of (12) can be evaluated. Similarly, evaluating the integrals in (11) leads to

The above two conditions furnish the last two necessary equations that are needed before the coefficient matrix can be stably inverted. In the next section, we study the invertibility of the coefficient matrix and the stability of the proposed numerical scheme.

3. Analysis of the Method

The above conditions provide () linearly independent equations that can be solved for the unknown coefficients , , , at each time increment. To prove the invertibility of the coefficient matrix we consider the case where . Other kernels can be treated in a similar way. It is appropriate to group the above conditions in the following order. It is possible to place (6) (which is evaluated at , or, ) in the first row. Next, one can place (7) (which is evaluated at ) in the second row. Following this pattern as well as evaluating two equations for every mid-point leads to equations that can be placed in the first rows of the matrix. The last two rows can be used to impose the boundary conditions given in (12) and (13). Following this manner, it is possible to group the above equations at every time interval in the form of

or

where , the unknown vector contains the constants of itegrations at every time interval, and is the known right hand side. In the above matrix equation, the constants on the right hand side are given by

and the parameters , and , are given by

It is now possible to denote and partition the above coefficient matrix according to

where the submatrices , , , , and are () matrices. The partitioned matrices are , , , and and are given by

The determinant of the coefficient matrix in (18) is given by

It is sufficient to show that the two determinants on the right hand side in the above equation are nonzero. The matrix can be inverted since it is an upper block matrix, and . In order to show that the second determinant is nonzero, one can denote and note that

It is also helpful to note that

Therefore, the second determinant on the right hand side of (20) is simplified to

where, for simplicity, , and

The above matrices are matrices, and one can obtain an expression for the determinant. It is possible to simplify the above determinant which leads to

By letting in the second summation, the above relation is further simplified to

Converting the rest of the terms to hyperbolic functions leads to

For a fixed , and the above determinant is equal to the difference between two positive monotonic functions of with different slopes. Therefore, there are infinite values of for which the determinant is not equal to zero.

Note that, for the case where the kernels in the boundary conditions, , are not constants, the matrix remains the same, that is, , and one needs to only check the in  (20).

In order to study the stability of the method, one notes that the right hand side of (14) includes the terms . Following (14), it follows that

The temperature fields can be written in terms of the coefficients . It follows that (14) can be written in terms of the coefficients at the previous time according to

where is given as before, and is given by

For stability, we require that all eigenvalues of the matrix be inside the unit circle. Figure 2 shows the magnitude of the largest eigenvalue as a function of the order of the approximation . As the mesh is refined, and , the magnitude of the largest eigenvalue increases. However, it stays within the unit circle, and the scheme is stable. The above stability test can also be applied to problems with variable kernel functions.

Figure 2: Norm of the largest eigenvalue of the matrix as a function of the order of the approximation .

4. Numerical Examples

Example 1. Consider the heat conduction given by (1) and (2) where

The exact solution for the above problem is given by . Choosing the values of and for the present method, Figure 3 presents the error at three locations within the domain. The method can produce numerical solution that compares well with the exact solution.

Figure 3: The error in the numerical solution at , , and for Example 1.

Example 2. Consider the heat conduction given by (1) and (2) where

This example was considered in [9, 10] to compare the standard finite difference approximations. The exact solution is given by for the initial condition given by . Using the values of , , and for the present method, Figure 4 compares the computed values of temperature at , , and to the exact values.
We next apply the present method to three nonlinear problems with nonlocal conditions.

Figure 4: The error in the numerical solution at , , and for Example 2, with .

Example 3. Consider the parabolic system given by

with

where , , and are given constants. The above problem involves a reaction between two reactants [16]. The nonlinear term appears in the boundary condition at . Applying the present method leads to the field variable given by

and formulating similar continuity conditions at mid-points leads to equations similar to (6) and (7). Equation (34) can also be treated similar to (9). The nonlinear boundary condition in (33) leads to

where . Note that in addition to , for , the scalar is also unknown. The nonlinear condition in (33) leads to quadratic nonlinear terms and in (36). A common approach is to set up an iteration at every time step according to

where is an initial guess, which is often the value of at the previous time interval. It is now possible to invert the coefficient matrix and solve for the , for and . One can use this value of for and repeat this process until the values converge. This is often the usual approach when one is faced with solving a nonlinear problem. We follow this approach in this example. Using the values of , , and from [16], and letting and , one can obtain the numerical results for the above problem. Figure 5 presents the numerical results using the present algorithm for and as a function of time. The results compare very well with the results reported in [16]. In particular, Table 1 compares the computed values of and for a few time instants. Note that the reported results in [16] are also approximate solutions and are reported after interpolation and retaining a few terms in the power series. We next consider a fully nonlinear problem and show the effectiveness of the present method. For a large class of nonlinear problems, the present method can provide a linear solution with no need for iterations at every time interval.

Table 1: Comparison of the computed values of and with the truncated results reported in [16] for = 100, = 0.01, and .

Time
Reference [16] truncated
Present
Reference [16] truncated
Present

= 0.4
0.6739
0.6877
0.0389
0.0419

= 1.0
0.8613
0.8788
0.0776
0.0798

= 2.4
0.9046
0.9047
0.0951
0.0950

Figure 5: Calculated values of and using present method for Example 3 from [16] for the case of , , . The comparison to the reported values is presented in Table 1.

Example 4. Consider the reaction-diffusion system given by [17]

where a numerical method based on finite difference, which involves iterations, was presented. Applying the present method it is possible to formulate a numerical method which does not need iterations and requires the inversion of only one matrix. Using the constant values of , and following the steps in (2) and (3) one can arrive at

Denoting and noting that , it is now possible to look for a solution in the form of . Considering various orders of leads to

Applying the same procedure to the boundary condition given in (38) leads to

and for various orders of

The above equations are linear, and applying the local formulation for various orders of approximation will require the inversion of only one matrix. The exact solution is given by

where

Dividing the domain into 400 intervals, that is, , and using , Figure 6 presents the absolute values of the error between the computed values and the exact solution. Note that the method requires the inversion of one matrix with dimension only. For this case and the first two terms in the expansion of are sufficient.

Figure 6: The error in the numerical solution at four points within the domain for Example 4.

Example 5. Consider the reaction-diffusion system given by [18]

For this system the quantity is conserved. The above model appears in the study of biological and chemical systems [18]. Our interest here is towards the numerical simulation of this model, and in particular the appearance of nonlinear terms. Following the approach presented here as well as writing the above equation in terms of local coordinates leads to

Denoting and noting that , it is possible to seek a solution in the form of , which leads to

Similar to the previous example, we only need to solve linear systems at each order of . In fact, they all require the inverse of the same matrix. As a result, the present method requires the inversion of only one matrix. Dividing the domain into 80 intervals, that is, , as well as using leads to . Depending on the initial condition, the above system can lead to blowup or a steady-state solution [18]. Consider the numerical solution of the above system when the initial condition is given by

Starting from the above initial condition, the steady state solution is a constant value of . Figure 7 shows the difference between the conserved value and its value at , that is, . The method requires the inversion of only one matrix.

Figure 7: The error in the conserved value as a function of time for Example 4.

5. Conclusion

In this paper we presented a local method for obtaining the solution for 1-D parabolic problems with nonlocal boundary conditions. The method is first order accurate in time. However, the solution is locally exact in the space dimension. The order of accuracy in time can easily be improved to second order. Five numerical examples were used to study the applicability of the method. In particular, three nonlinear problems were studied to show the strength of the method.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

J. R. Cannon and J. van der Hoek, “Diffusion subject to the specification of mass,” Journal of Mathematical Analysis and Applications, vol. 115, no. 2, pp. 517–529, 1986.

A. Bouziani, “On the quasi static flexure of a thermoelastic rod,” Communications in Applied Analysis, vol. 6, no. 4, pp. 549–568, 2002.

W. A. Day, “Extensions of a property of the heat equation to linear thermoelasticity and other theories,” Quarterly of Applied Mathematics, vol. 40, no. 3, pp. 319–330, 1982.

Y. S. Choi and K. Y. Chan, “A parabolic equation with nonlocal boundary conditions arising from electrochemistry,” Nonlinear Analysis. Theory, Methods & Applications, vol. 18, no. 4, pp. 317–331, 1992.

V. Capasso and K. Kunisch, “A reaction-diffusion system arising in modelling man-environment diseases,” Quarterly of Applied Mathematics, vol. 46, no. 3, pp. 431–449, 1988.

J. Furter and M. Grinfeld, “Local vs. non-local interactions in population dynamics,” Journal of Mathematical Biology, vol. 27, no. 1, pp. 65–80, 1989.

T.-T. Li, “A class of nonlocal boundary value problems for partial differential equations and its applications in numerical analysis,” vol. 28, pp. 49–62.

A. Štikonas, “A survey on stationary problems, Green's functions and spectrum of Sturm-Liouville problem with nonlocal boundary conditions,” Nonlinear Analysis: Modelling and Control, vol. 19, no. 3, pp. 301–334, 2014.

G. Ekolin, “Finite difference methods for a nonlocal boundary value problem for the heat equation,” BIT: Numerical Mathematics, vol. 31, no. 2, pp. 245–261, 1991.

Y. Liu, “Numerical solution of the heat equation with nonlocal boundary conditions,” Journal of Computational and Applied Mathematics, vol. 110, no. 1, pp. 115–127, 1999.

L. Mu and H. Du, “The solution of a parabolic differential equation with non-local boundary conditions in the reproducing kernel space,” Applied Mathematics and Computation, vol. 202, no. 2, pp. 708–714, 2008.

L. Bougoffa and R. Rach, “Solving nonlocal initial-boundary value problems for linear and nonlinear parabolic and hyperbolic PDEs by the Adomian expansion,” Applied Mathematics and Computation, vol. 225, pp. 50–61, 2013.

M. Tadi, “A fixed-grid local method for 1-D Stefan problems,” Applied Mathematics and Computation, vol. 219, no. 4, pp. 2331–2341, 2012.

J. Zhou and C. Mu, “Blowup for a degenerate and singular parabolic equation with non-local source and absorption,” Glasgow Mathematical Journal, vol. 52, no. 2, pp. 209–225, 2010.

Y. Chen, Q. Liu, and C. Xie, “Blow-up for degenerate parabolic equations with nonlocal source,” Proceedings of the American Mathematical Society, vol. 132, no. 1, pp. 135–145, 2004.

J. A. Dixon, “A nonlinear weakly singular Volterra integro-differential equation arising from a reaction-diffusion study in a small cell,” Journal of Computational and Applied Mathematics, vol. 18, no. 3, pp. 289–305, 1987.

C. V. Pao, “Numerical solutions of reaction-diffusion equations with nonlocal boundary conditions,” Journal of Computational and Applied Mathematics, vol. 136, no. 1-2, pp. 227–243, 2001.

C. Budd, B. Dold, and A. Stuart, “Blowup in a partial differential equation with conserved first integral,” SIAM Journal on Applied Mathematics, vol. 53, no. 3, pp. 718–742, 1993.