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

Authors: M. Tadi; Miloje Radenkovic
Journal: International Journal of Computational Mathematics
 (2014)

---

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


---

## Body

## 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].Section2 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(1)ut=uxx+qt,x,x∈0,1,u(0,x)=f(x),ut,0=∫01‍ϕxut,xdx,u(t,1)=∫01‍ψ(x)u(t,x)dx,
where the term q(t,x) is the heat generation, f(x) is the initial condition, and the boundary conditions are of nonlocal type. Assume that the domain is divided into m equal intervals. Consider a finite difference mesh given in Figure 1.Figure 1
Finite difference mesh showing the local coordinate.For a typical nodek, consider an implicit expansion of the equation at time interval (n+1) given by
(2)u-uknΔt=uxx+q,
where the value of ukn is known. Now, consider a transformation to the local coordinate, η, according to
(3)η=(x-xk)Δx.
In terms of the local coordinate, (2) simplifies to
(4)uηη-λ2u=-λ2ukn-Δx2qk,η∈-1,1,with,λ2=Δx2Δt.
The variable qk denotes the heat generation function q(t,x) evaluated at time (Δt(n+1)) and location k. (In general, it is a function of space η; however, for a first-order accurate method, the value at x=(k-1)Δx can be used. Using Taylor series expansion, it is possible to obtain a more accurate description in terms of η; that is, ∑l=0∞‍alηl, for which exact solution can be obtained.) Treating qk as alocal constant inη function allows us to write down the exact solution of u(t,η) around every node given by
(5)ut,η=ckeλη+dke-λη+ξkn+1,ξkn+1=ukn+Δtqkn+1,
where now ck and dk are constants of integration. For m intervals, there are (m+1) nodes. Therefore, there are M=m+1 equations similar to (5), which describe the temperature profile; that is, k=1,2,…,M. It follows that there are 2M constants, that is, ck,dk (k=1,2,…,M) that are unknown and should be computed at each time interval. It is possible to introduce 2M linearly independent conditions. These conditions are not unique. However, a linearly independent set of conditions is given by the following.(i)
The temperature distribution att(n+1)=(n+1)Δt should be continuous. We can impose this condition at mid-points between nodes. It leads to
(6)cke1/2λ+dke-1/2λ+ξkn+1=ck+1e-1/2λ+dk+1e1/2λ+ξk+1n+1,hhhhhhhhhhhhhhfor,k=1,…,m.(ii)
The flux att(n+1) should be continuous. In the present case, it follows that the temperature at t(n+1) should have a continuous first derivative. We can also impose this condition at mid-points between nodes. It leads to
(7)cke1/2λ-dke-1/2λ=ck+1e-1/2λ-dk+1e1/2λ,for,k=1,…,m.The above conditions furnish (2m) linearly independent equations. The additional two conditions needed are given by the boundary conditions. The solution must satisfy the boundary conditions at x=0 and x=1. These conditions are nonlocal and can be imposed according to
(8)u(tn+1,0)=c1+d1+ξ1n+1=∫01‍ϕ(x)u(tn+1,x)dx(9)u(tn+1,1)=cM+dM+ξMn+1=∫01‍ψ(x)u(tn+1,x)dx.
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 [0,Δx/2]. For the nodes k=2,3,…,m, the domains can be chosen to be [(k-1)Δx-Δx/2,(k-1)Δx+Δx/2]. For the last node, the domain can be chosen as [1-Δx/2,1]. In terms of the local coordinates η, for the interior nodes, the individual domains of integrations are all [-1/2,1/2]. Similarly, the domain of integration is [0,1/2] for the first node and [-1/2,0] for the last node. Using these domains, the integral condition at x=0 is given by
(10)c1+d1+ξ1n+1=∫01‍ϕ(x)u(tn+1,x)dx=Δx∫01/2‍ϕηc1eλη+d1e-λη+ξ1n+1dη+∑k=2m‍Δx∫-1/21/2‍ϕηckeλη+dke-λη+ξkn+1dη+Δx∫-1/20‍ϕηcMeλη+dMe-λη+ξMn+1dη,
and for x=1, it is given by
(11)cM+dM+ξMn+1=∫01‍ψ(x)u(tn+1,x)dx=Δx∫01/2‍ψηc1eλη+d1e-λη+ξ1n+1dη+∑k=2m‍Δx∫-1/21/2‍ψηckeλη+dke-λη+ξkn+1dη+Δx∫-1/20‍ψηcMeλη+dMe-λη+ξMn+1dη.
The functions ϕ(x) and ψ(x) are known, and the above integrals can be evaluated in closed forms. After evaluating the integrals, (10) leads to
(12)1-Δx∫01/2‍ϕηeληdηc1+1-Δx∫01/2‍ϕ(η)e-ληdηd1-∑k=2m‍Δx∫-1/21/2‍ϕηeληdηck-∑k=2m‍Δx∫-1/21/2‍ϕ(η)e-ληdηdk-Δx∫-1/20‍ϕηeληdηcM-Δx∫-1/20‍ϕηeληdηdM=-ξ1n+1+∫01‍ϕ(x)ξn+1dx.
Note that ξn+1 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
(13)-Δx∫01/2‍ψ(η)eληdηc1-Δx∫01/2‍ψ(η)e-ληdηd1-∑k=2m‍Δx∫-1/21/2‍ψ(η)eληdηck-∑k=2m‍Δx∫-1/21/2‍ψ(η)e-ληdηdk+Δx1-∫-1/20‍ψ(η)eληdηcM+Δx1-∫-1/20‍ψ(η)eληdηdM=-ξ1n+1+∫01‍ψ(x)ξn+1dx.
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 (2M) linearly independent equations that can be solved for the unknown coefficients ck, dk, k=1,2,…,M, at each time increment. To prove the invertibility of the coefficient matrix we consider the case where ϕ(x)=ψ(x)=1. 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 η=1/2, or, x=Δx/2) in the first row. Next, one can place (7) (which is evaluated at η=1/2) in the second row. Following this pattern as well as evaluating two equations for every mid-point leads to (2m) equations that can be placed in the first (2m) 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 tn+1 in the form of
(14)Avn+1=gn+1,
or
(15)eθe-θeθ-e-θ-e-θ-eθ-e-θeθ0000…0000eθe-θeθ-e-θ-e-θ-eθ-e-θeθ…………………………0000eθe-θeθ-e-θ-e-θ-eθ-e-θeθμ10ν10μ11ν11ϵϵϵϵ…μM0νM0μM1νM1×c1d1c2d2c3d3⋮cmdmcMdM=g10g20g30⋮gm0σδ,
where θ=λ/2, the unknown vector vn+1 contains the constants of itegrations at every time interval, and gn+1 is the known right hand side. In the above matrix equation, the constants on the right hand side are given by
(16)gk=ξk+1n+1-ξkn+1,fork=1,2,…,mσ=-1+Δx2ξ1n+1+Δx∑k=2m‍ξkn+1+Δx2ξMn+1δ=Δx2ξ1n+1+Δx∑k=2m‍ξkn+1+Δx2-1+Δx2ξMn+1,
and the parameters μ10, μ11 and νM0, νM1 are given by
(17)μ10=1-Δxλeθ-1,ν10=1+Δxλ(e-θ-1),μM0=-Δxλ(1-e-θ),νM0=Δxλ(1-eθ),ϵ=Δxλ(e-θ-eθ),μ11=-Δxλeθ-1,ν11=Δxλ(e-θ-1),μM1=1-Δxλ1-e-θ,νM1=1+Δxλ(1-eθ).
It is now possible to denote and partition the above coefficient matrix according to
(18)DC0000DC00⋯⋯⋯⋯⋯000DCE1F…FEm,or,Γ11Γ12Γ21Γ22,
where the submatrices D, C, E1, Em, and F are (2×2) matrices. The partitioned matrices are Γ11∈R(M-2)×(M-2), Γ22∈R(2×2), Γ12∈R(M-2)×2, and Γ22∈R(2×2) and are given by
(19)Γ11=DC000DC0⋯⋯⋯⋯000D,Γ12=00⋯C,Γ21=E1F…F,Γ22=Em.
The determinant of the coefficient matrix in (18) is given by
(20)det⁡(A)=det⁡(Γ11)︸Θ1det⁡(Γ22-Γ21Γ11-1Γ12)︸Θ2.
It is sufficient to show that the two determinants on the right hand side in the above equation are nonzero. The matrix Γ11 can be inverted since it is an upper block matrix, and det⁡(D)=-2≠0. In order to show that the second determinant is nonzero, one can denote G=D-1C and note that
(21)Γ11-1=D-1-[GD-1][G2D-1]⋯0D-1-[GD-1][G2D-1]⋯⋯D-1-[GD-1]000D-1.
It is also helpful to note that
(22)G=-H=-e-λ00eλ,Gl=-1lHl=-1le-lλ00elλ.
Therefore, the second determinant on the right hand side of (20) is simplified to
(23)Θ2=det⁡(Γ22-Γ21Γ11-1Γ12)=det⁡Em+E1Hm+F∑j=1m-1‍Hj,Θ2=det⁡μM0+d1μ10νM0+ν11d2μM1+d1μ11νM1+ν11d2+ϵc1c2c1c2,
where, for simplicity, ϵ=-Δx/λ(eθ-e-θ)=-2Δx/λsinh⁡(θ), and
(24)d1=e-mλ,d2=emλ,d1d2=1,c1=∑l=1m-1‍e-lλ,c2=∑l=1m-1‍elλ.
The above matrices are (2×2) matrices, and one can obtain an expression for the determinant. It is possible to simplify the above determinant which leads to
(25)Θ2=2Δt(eθ+e-θ-2)+4Δtsinh⁡(θ)∑l=1m-1‍sinh⁡(lλ)+4Δtsinh⁡(θ)∑l=1m-1‍sinh⁡((m-l)λ)+1-2Δteθ-1e-mλ-1+2Δte-θ-1emλ.
By letting k=m-l in the second summation, the above relation is further simplified to
(26)Θ2=2Δt(eθ+e-θ-2)+8Δtsinh⁡(θ)∑l=1m-1‍sinh⁡(lλ)+1-2Δteθ-1e-mλ-1+2Δte-θ-1emλ.
Converting the rest of the terms to hyperbolic functions leads to
(27)Θ2=4Δtcosh⁡λ2+8Δtsinh⁡λ2∑l=1m-1‍hh×∑l=1m-1‍sinh⁡⁡lλ+4Δtcosh⁡⁡mλ-4Δt+2sinh⁡⁡mλ+4Δtcosh⁡⁡m+12λ.
For a fixed Δt, λ=Δx/Δt and the above determinant is equal to the difference between two positive monotonic functions of Δx with different slopes. Therefore, there are infinite values of Δx for which the determinant is not equal to zero.Note that, for the case where the kernels in the boundary conditions,ϕ(x),ψ(x), are not constants, the matrix Γ11 remains the same, that is, (det⁡(Γ11)≠0), and one needs to only check the Θ2 in  (20).In order to study the stability of the method, one notes that the right hand side of (14) includes the terms gk. Following (14), it follows that
(28)gk=ξk+1n+1-ξkn+1=uk+1n-ukn+Δtqk+1n+1-qk+1n.
The temperature fields ukn can be written in terms of the coefficients ukn=ckn+dkn+ξkn. It follows that (14) can be written in terms of the coefficients at the previous time according to
(29)Avn+1=Bvn+h,
where A is given as before, and B is given by
(30)Avn+1=-1-10011000000…0000-1-100100…………………………0000-1-100110000000000…0000×c1d1c2d2c3d3⋮cmdmcMdM+g10g20g30⋮gm0σδ.
For stability, we require that all eigenvalues of the matrix A-1B be inside the unit circle. Figure 2 shows the magnitude of the largest eigenvalue as a function of the order of the approximation m. As the mesh is refined, and Δx→0, 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 matrixA-1B as a function of the order of the approximation m.


## 4. Numerical Examples

Example 1.
Consider the heat conduction given by (1) and (2) where
(31)q(t,x)=(1+x)(2e-t-1-te-t),ϕ(x)=1.2x,ψ(x)=4/3.
The exact solution for the above problem is given by u(t,x)=(1+x)(1-e-t-t+te-t). Choosing the values of Δx=0.01 and Δt=0.00001 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 atx=0.3, x=0.5, and x=0.7 for Example 1.Example 2.
Consider the heat conduction given by (1) and (2) where
(32)q(t,x)=-e-tx(x-1)+δ6(1+δ)+2,ϕ(x)=-δ,ψ=-δ.
This example was considered in [9, 10] to compare the standard finite difference approximations. The exact solution is given by e-tx(x-1)+δ/6(1+δ) for the initial condition given by f(x)=x(x-1)+δ/6(1+δ). Using the values of δ=0.0144, Δt=0.0025, and Δx=0.0025 for the present method, Figure 4 compares the computed values of temperature at x=0.32, x=0.5, and x=0.64 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 atx=0.32, x=0.5, and x=0.64 for Example 2, with δ=0.0144.Example 3.
Consider the parabolic system given by(33)ut=uxx,0<x<1,t>0,ux(t,0)=0,t>0ux(t,1)=Em1+LLθ(t)-(1-θ(t))u(t,1)u(0,x)=1,0<x<1,θ(0)=0,
with
(34)mθ(t)+∫01‍u(t,x)dx=1,t>0,
where E, L, and m are given constants. The above problem involves a reaction between two reactants [16]. The nonlinear term appears in the boundary condition at x=1. Applying the present method leads to the field variable given by
(35)u(t,η)=ckeλη+dke-λη+ukn,
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
(36)(λ+β)cM+(-λ+β)dM-β(L-uMn)θ-βθcM-βθdM=-βuMn,
where β=EmΔx/1+L. Note that in addition to ck,dk, for k=1,…,M, the scalar θ(t) is also unknown. The nonlinear condition in (33) leads to quadratic nonlinear terms θcM and θdM in (36). A common approach is to set up an iteration at every time step according to
(37)(λ+β-βθ^)cM+(-λ+β-βθ^)dM-β(L-uMn)θ=-βuMn,
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 ck,dk, for k=1,…,M 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 E=100, m=1, and L=0.01 from [16], and letting Δx=1/50 and Δt=0.0005, one can obtain the numerical results for the above problem. Figure 5 presents the numerical results using the present algorithm for θ(t) and u(t,1) 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 θ(t) and u(t,1) 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θ(t) and u(t,1) with the truncated results reported in [16] for E = 100, L = 0.01, and m=1.0.
Time
Reference [16] truncated θ
Present
Reference [16] truncated u(t,1)
Present
t = 0.4
0.6739
0.6877
0.0389
0.0419
t = 1.0
0.8613
0.8788
0.0776
0.0798
t = 2.4
0.9046
0.9047
0.0951
0.0950Figure 5
Calculated values ofθ(t) and u(t,1) using present method for Example 3 from [16] for the case of E=100, L=0.01, m=1. The comparison to the reported values is presented in Table 1.Example 4.
Consider the reaction-diffusion system given by [17]
(38)ut=Duxx-au+u1-u+qt,x,0<x<1,t>0,ut,0=0,t>0ux(t,1)+u(t,1)=12∫01‍u(t,x)dx+g(t),u(0,x)=x(2-x)+sinπx2,
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 D=0.1, a=5 and following the steps in (2) and (3) one can arrive at
(39)uηη-Δx240+10Δtu-10Δx2u2=-10u^Δt+10q(t,xi)Δx2.
Denoting λ2=Δx240+10/Δt and noting that ϵ=10Δx2≪1, it is now possible to look for a solution in the form of u=u0+ϵu1+ϵ2u2+⋯. Considering various orders of ϵ leads to
(40)ϵ0,d2u0dη2-λ2u0=-10u^Δt+10q(t,xi)Δx2=ϕkϵ1,d2u1dη2-λ2u1=u02ϵ2,d2u2dη2-λ2u2=2u0u1⋮
Applying the same procedure to the boundary condition given in (38) leads to
(41)uη+Δxu=12Δx∫01‍udx+Δxg(t)
and for various orders of ϵ(42)ϵ0,du0dηx=1+Δxu0x=1-12Δx∫01‍u0dx=Δxg(t)ϵ1,du1dηx=1+Δxu1x=1-12Δx∫01‍u1dx=0ϵ2,du2dηx=1+Δxu1x=1-12Δx∫01‍u2dx=0.
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
(43)u*t,x=x2-x+0.1e-αtsinπx2,withα=5+π240,
where
(44)g(t)=ux*x=1+u*x=1-12∫01‍u*(t,x)dx.
Dividing the domain into 400 intervals, that is, Δx=1/400, and using Δt=0.0005, 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 (401×401) only. For this case ϵ=0.000625 and the first two terms in the expansion of u(η) are sufficient.Figure 6
The error in the numerical solution at four points within the domain[0,1] for Example 4.Example 5.
Consider the reaction-diffusion system given by [18]
(45)ut=uxx+u2-∫01‍u2dx,0<x<1,t>0,u0,x=ψx,ux(t,0)=ux(t,1)=0.
For this system the quantity ∫01‍u(t,x)dx 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
(46)uηη-Δx2Δtu+Δx2u2-Δx2∫01‍u2dx=-Δx2Δtukn.
Denoting λ2=Δx2/Δt and noting that ϵ=Δx2≪1, it is possible to seek a solution in the form of u=u0+ϵu1+ϵ2u2+⋯, which leads to
(47)ϵ0,d2u0dη2-λ2u0=-λ2uknϵ1,d2u1dη2-λ2u1=-u02+∫01‍u02dxϵ2,d2u2dη2-λ2u2=-2u0u1+2∫01‍u0u1dx⋮
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, Δx=1/100, as well as using Δt=0.000125 leads to ϵ=0.0001. 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
(48)ψ(x)=e-(x-0.5)4/0.004.
Starting from the above initial condition, the steady state solution is a constant value of u=0.4559=∫01‍ψ(x)dx. Figure 7 shows the difference between the conserved value ∫01‍u(t,x)dx and its value at t=0, that is, ∫0ψ‍(x)dx. The method requires the inversion of only one matrix.Figure 7
The error in the conserved value as a function of time for Example4.


## 5. Conclusion

In this paper we presented alocal 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.



---

*Source: 923693-2014-11-20.xml*