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.
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.
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.
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.
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.
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.
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.
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.