# Highly Efficacious Sixth-Order Compact Approach with Nonclassical Boundary Specifications for the Heat Equation

Authors: Nouria Arar; Leila Ait Kaki; Abdellatif Ben Makhlouf
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper suggests an accurate numerical method based on a sixth-order compact difference scheme and explicit fourth-order Runge–Kutta approach for the heat equation with nonclassical boundary conditions (NCBC). According to this approach, the partial differential equation which represents the heat equation is transformed into several ordinary differential equations. The system of ordinary differential equations that are dependent on time is then solved using a fourth-order Runge–Kutta method. This study deals with four test problems in order to provide evidence for the accuracy of the employed method. After that, a comparison is done between numerical solutions obtained by the proposed method and the analytical solutions as well as the numerical solutions available in the literature. The proposed technique yields more accurate results than the existing numerical methods.


---

## Body

## 1. Introduction

Many mathematical patterns of real-world problems engender partial differential equations (PDEs) of initial and boundary conditions. Each of these issues can be represented in different ways by a hyperbolic, elliptical, or parabolic partial differential equation. They can be homogeneous or inhomogeneous, in one, two, or three dimensions, either with local or nonlocal boundary conditions besides the initial conditions. Different researchers have dealt with partial differential equations with nonlocal boundary conditions. Most of the studies were concerned with the second-order parabolic equation, mainly the heat equation. The research of the numerical solution for the nonclassical boundary value problems has recently been an important topic of research in various fields of science and engineering, such as thermoelasticity [1]. Some issues appearing in thermodynamics [2], heat conduction [3–10], and plasma physics [11] can be reduced to the nonclassical problems with the integral condition.In this paper, the following one-dimensional problem is considered:(1)∂w˜∂τ=∂2w˜∂ξ2+fξ,τ,0<ξ<1,0<τ≤T,where fξ,τ is a sufficiently differentiable function in time and space, T∈R+, and subject to the initial condition(2)w˜ξ,0=φξ,0≤ξ≤1,and the NCBC(3)w˜0,τ=∫01f1ξw˜ξ,τdξ+g1τ,0≤τ≤T,w˜1,τ=∫01f2ξw˜ξ,τdξ+g2τ,0≤τ≤T,where f1, f2, g1, g2, and φ are known functions.Many authors examined this type of parabolic initial-boundary value problems in one dimension involving nonclassical boundary conditions; some examples of these problems can be mentioned as follows:(i)
Without doing any numerical experiments, Cannon and van der Hoek [12] developed an implicit finite difference approach for the problem given by (1), when fξ,τ=0, with boundary condition(4)w˜1,τ=gτ,0≤τ≤1,and subject to the energy specification(5)∫0bw˜ξ,τdξ=Ξτ,0<τ≤T,0<b<T.Deckert et al. demonstrated the existence and uniqueness of the solutions using the Laplace transform in [13].(ii)
Makarov and Kulyev in [14] considered the problem (1), when the term source depends on the solution, i.e., fξ,τ is replaced by fξ,τ,w˜ with homogeneous initial condition and the boundary condition(6)w˜0,τ=0,0<τ<T,as well as the NCBC(7)∂w˜∂ξ0,τ=∂w˜∂ξ1,τ,0<τ<T.(iii)
Shingmin and Yanping in [15] chose the Crank-Nicolson method and Simpson’s numerical integration formula to solve the problem (1) with initial condition (2), and nonlocal boundary conditions(8)w˜0,τ=∫01ϕξ,τw˜ξ,τdξ,0<τ≤T,w˜1,τ=∫01ψξ,τw˜ξ,τdξ,0<τ≤T,Information about the existence and uniqueness of the solution is provided in [16].(iv)
Ekolin in [17] expanded and assessed three possible finite difference methods to solve the problem (1), with an initial condition (2) and integral boundary conditions (3). The schemes advanced in [17] are based on the forward and backward Euler techniques and the Crank-Nicolson method. The same problem is also considered by Ang in [18] where the author proposed another method of the solution not based on finite difference approximations.(v)
Cannon and Matheson in [19] developed a numerical scheme for the problem (1) with initial condition (2), boundary condition (4), and the NCBC (5), when the upper bound of the integral in (5) depends on τ, i.e., 0<b=bτ<T. The existence, uniqueness, and computation of the solution can be consulted in [20].(1) Liu [21] advocated the θ-schema to solve numerically the problem (1) with an initial condition (2) and nonlocal boundary conditions (3).(2) Bastani and Salkuyeh [22] treated the problem(9)∂w˜∂τ−α∂2w˜∂ξ2=fξ,τ,a<ξ<b,0<τ≤T,with nonlocal boundary condition (5), such that ξ vary from a to b, the initial condition (2), and the boundary condition(10)w˜0,τ=gτ,0<τ≤T.(vi)
By using a spectral collocation method with preconditioning and cubic interpolation for approximating the nonlocal boundary condition, El-Gamel et al. in [23] proposed a highly efficient algorithm for solving the parabolic problem (1) with initial condition (2) and nonlocal boundary conditions(11)∑p=01αpτ∂p∂ξpw˜0,τ+∫01φ0ξ,τw˜x,τdξ=ψ0t,0<τ≤T,∑p=01βpτ∂p∂ξpw˜1,τ+∫01φ1ξ,τw˜ξ,tdξ=ψ1τ,0<τ≤T,in which they employ orthogonal Chelyshkov polynomials as the basis.The remaining sections are arranged as follows: the high accuracy compact finite difference formula for the second-order derivative is briefly introduced in Section2 of this document. Section 3 describes how to solve the heat equation with integral boundary conditions using the compact finite difference method and explicit fourth-order Runge-Kutta approach. Numerical experiments are presented in Section 4 to corroborate our findings. Section 5 draws an appropriate conclusion.To develop our proposed technique, we shall in what follows divide the spatial interval0,1 into ℵ+1 subintervals having length h=1/ℵ+1, and the time interval is divided into subintervals having length Δτ. The grid points are detailed by(12)ξi=ih,i=0,ℵ+1¯,τj=jΔτ,j=0,ℵτ¯.ℵ and ℵτ are integers. Here, ℵ is chosen as an odd number because of the use of Simpson’s one-third rule formula.It is noteworthy that we useϑξ,τ to denote the numerical solution for the exact solution w˜ξ,τ of problem (1), in the following sections.


## 2. Sixth Compact Scheme for Second-Order Derivative

The most well-known numerical method for approximating differential equation solutions is the finite difference method (FDM). Higher-order FDMs play a vital role in the numerical resolution of various PDEs; however, when the schemes are explicit, they require the use of a large number of stencils leading to the creation of large matrix bandwidths, and so complicating the numerical processing. Implicit by construction, compact FD schemes like those found [24–26] are recognised to have better spectral resolution and lead to high-order accuracy by using additional derivative information. Furthermore, the implicit high-order compact schemes are very flexible to use because they are mainly composed of fixed-difference equations, relating the nodal values of the original function and its derivatives. Since the topic of interest can be changed, there is no need to exercise extra caution because these different equations are independent of the primary equation.Since 1970, compact finite difference methods have been considerably used to solve a large number of differential equations such as the advection-diffusion equation [27], Bateman–Burgers equation [28], Black–Scholes equation [29], diffusion equation [30], integro-differential equation [31], singular boundary problems [32], and wave equation [33].In this section, we give a development of the sixth-order compact finite difference formula for the second-order derivative. Lele in [26] defined compact schemes for approximating second-order derivatives at interior nodes as(13)lϑi−1″+ϑi″+lϑi+1″=pϑi+2−2ϑi+ϑi−24h2+qϑi+1−2ϑi+ϑi−1h2,i=2,ℵ−1.¯Forl=2/11, the scheme is sixth-order accurate, and in this case, q=12/11, p=3/11, and the truncation error is −4/6!11l−2h4ϑ6.By simplifying the relation (13), we obtain(14)lϑi−1″+ϑi″+lϑi+1″=m1ϑi−2+m2ϑi−1−m3ϑi+m2ϑi+1+m1ϑi+2,where(15)m1=p4h2,m2=qh2,m3=p+4q2h2.In order to get the same accuracy and dominant error term, special formulas are needed fori=1 and i=ℵ. When i=1, we use(16)ϑ1″+l1ϑ2″=1h2s0ϑ0+s1ϑ1+s2ϑ2+s3ϑ3+s4ϑ4+s5ϑ5+s6ϑ6.Wheni=ℵ, we use(17)ϑℵ″+l1ϑℵ−1″=1h2s0ϑℵ+1+s1ϑℵ+s2ϑℵ−1+s3ϑℵ−2+s4ϑℵ−3+s5ϑℵ−4+s6ϑℵ−5.where the coefficients can be found by matching the Taylor series expansions up to the order of h7, which gives us a linear system in which the solution is given by(18)l1=112,s0=131360,s2=−574,s3=15718,s4=−98,s5=320,s6=−190.All formulas areOh6. More details are available in [24–26]. When unknowns ϑ0 and ϑℵ+1 are managed by the boundary conditions of the considered problem, we obtain a system of ℵ equations with ℵ unknowns ϑ1, …, ϑℵ.


## 3. Handling Numerical Scheme for the Heat Equation

The development of a sixth-order compact finite difference formula is presented in this section in order to obtain a numerical solution to the nonhomogeneous diffusion equation with NCBC. Sixth-order finite difference approximations are used to approximate the spatial derivative. Simpson’s one-third rule formula is used to meet the difficulty of the integral boundary conditions while also eliminating additional variables to obtain a system ofℵ equations with ℵ variables.The combination of (1) with formulas (14), (16), and (17) at time level τ=τj constructs a system of ℵ linear equations with ℵ+2 unknowns ϑ0, ϑ1, …, ϑℵ+1. Simpson’s one-third rule formula is used to approximate the integrals in (3) and thus eliminate ϑ0 and ϑℵ+1, to obtain a system of ℵ ordinary differential linear equations in ℵ unknowns ϑ1, …, ϑℵ.(19)ϑ0τ=h3f1ξ0ϑ0τ+f1ξℵ+1ϑℵ+1τ+2∑ι=1ℵ+1/2−1f1ξ2ιϑξ2ι,τ+4∑ι=1ℵ+1/2f1ξ2ι−1ϑξ2ι−1,τ+g1τ,ϑℵ+1τ=h3f2ξ0ϑ0τ+f2ξℵ+1ϑℵ+1τ+2∑ι=1ℵ+1/2−1f2ξ2ιϑξ2ι,τ+4∑ι=1ℵ+1/2f2ξ2ι−1ϑξ2ι−1,τ+g2τ.So, we have(20)ϑ0τ=h3D2∑ι=1ℵ+1/2−1η2f1ξ2ι+μ1f2ξ2ιϑξ2ι,τ+4∑ι=1ℵ+1/2η2f1ξ2ι−1+μ1f2ξ2ι−1ϑξ2ι−1,τ+η2g1τ+μ1g2τD,(21)ϑℵ+1τ=h3D2∑ι=1ℵ+1/2−1μ2f1ξ2ι+η1f2ξ2ιϑξ2ι,τ+4∑ι=1ℵ+1/2μ2f1ξ2ι−1+η1f2ξ2ι−1ϑξ2ι−1,τ+μ2g1τ+η1g2τD,where(22)η1=1−h3f1ξ0,η2=1−h3f2ξℵ+1,μ1=h3f1ξℵ+1,μ2=h3f2ξ0,D=η1⋅η2−μ1⋅μ2.


### 3.1. Solving Heat Equation by Compact Finite Difference Method

Combining (1)with formulas (14), (16) and (17), and (20) and (21) generates ℵ×ℵ system of linear ordinary differential equations which can be written in vector matrix form as(23)dVdτ=A−1MVτ+Lτ+Fτ,V0=Φξ,where(i)
V=ϑ1,ϑ2,…,ϑℵT is the approximate solution of heat equation (1)(ii)
Φξ=φξ1,…,φξℵT is the primary dataLet us, respectively, denoteg1τ and g2τ by g1 and g2, and then, the vectors L and F are given by(24)Lτ=1Ds0η2g1+μ1g2h2,m1η2g1+μ1g2,0,…,0,m1μ2g1+η1g2,s0μ2g1+η1g2h2T.And(25)Fτ=fξ1,τ,fξ2,τ,…,fξℵ,τT.A and M are matrices resulting from the sixth-order compact difference scheme defined by(26)A=1l1000…0l1l00…00l1l⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱l1l00…00l1l0…000l11.(27)M=M1,1M1,2M1,3M1,4M1,5M1,6…M1,ℵM2,1M2,2M2,3M2,4M2,5M2,6…M2,ℵm1m2−m3m2m10…00m1m2−m3m2m1⋱⋮⋮⋱⋱⋱⋱⋱⋱00…0m1m2−m3m2m1Mℵ−1,1…Mℵ−1,ℵ−5Mℵ−1,ℵ−4Mℵ−1,ℵ−3Mℵ−1,ℵ−2Mℵ−1,ℵ−1Mℵ−1,ℵMℵ,1…Mℵ,ℵ−5Mℵ,ℵ−4Mℵ,ℵ−3Mℵ,ℵ−2Mℵ,ℵ−1Mℵ,ℵ,where for i=1,3¯(28)M1,2i−1=s2i−1+s04h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M1,2i=s2i+s02h3Dη2f1ξ2i+μ1f2ξ2i.Fori=4,ℵ+1/2−1¯(29)M1,2i−1=s04h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M1,2i=s02h3Dη2f1ξ2i+μ1f2ξ2i,M1,ℵ=s04h3Dη2f1ξℵ+μ1f2ξℵ.The last line for the matrixM is given by: for i=0,2¯(30)Mℵ,ℵ−2i=s2i+1+s04h3Dμ2f1ξℵ−2i+η1f2ξℵ−2i,Mℵ,ℵ−2i−1=s2i+2+s02h3Dμ2f1ξℵ−2i−1+η1f2ξℵ−2i−1,Mℵ,ℵ−6=s04h3Dμ2f1ξℵ−6+η1f2ξℵ−6.Fori=1,ℵ−5/2−1¯,(31)Mℵ,2i−1=s04h3Dμ2f1ξ2i−1+η1f2ξ2i−1,Mℵ,2i=s02h3Dμ2f1ξ2i+η1f2ξ2i.The second line of the matrixM is given by(32)M2,i=m2+m14h3Dη2f1ξi+μ1f2ξi,fori=1,3,M2,2=−m3+m12h3Dη2f1ξ2+μ1f2ξ2,M2,4=m1+m12h3Dη2f1ξ4+μ1f2ξ4.Fori=3,ℵ+1/2−1¯,(33)M2,2i−1=m14h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M2,2i=m12h3Dη2f1ξ2i+μ1f2ξ2i,M2,ℵ=m14h3Dη2f1ξℵ+μ1f2ξℵ.The lineℵ−1 for the matrix M is given by(34)Mℵ−1,i=m2+m14h3Dμ2f1ξi+η1f2ξi,fori=ℵ,ℵ−2,Mℵ−1,ℵ−1=−m3+m12h3Dμ2f1ξℵ−1+η1f2ξℵ−1,Mℵ−1,ℵ−3=m1+m12h3Dμ2f1ξℵ−3+η1f2ξℵ−3,Mℵ−1,ℵ−4=m14h3Dμ2f1ξℵ−4+η1f2ξℵ−4.Fori=1,ℵ−3/2−1¯,(35)Mℵ−1,2i−1=m14h3Dμ2f1ξ2i−1+η1f2ξ2i−1,Mℵ−1,2i=m12h3Dμ2f1ξ2i+η1f2ξ2i.


### 3.2. Resolution of Ordinary Differential System (23)

In order to solve the obtained initial-value problem for a system ofN linear ordinary differential equations (23), we consider the single-step difference schemes for calculating the numerical values of the solution. Let uniform grid space for interval 0,T given by τ0=0 and τj=jΔτ for all j=0,1,…,ℵτ, where Δτ>0 is constant. This difference scheme is a gradual method; we start from the given approximate values Vτj−1 and proceed stepwise, calculating approximate values Vτj at the mesh points. The system (23) is solved using the explicit fourth-order Runge–Kutta scheme which is formulated by(36)k1=Gτj−1,Vτj−1,k2=Gτj−1+12Δτ,Vτj−1+12Δτk1,k3=Gτj−1+12Δτ,Vτj−1+12Δτk2,k4=Gτj−1+Δτ,Vτj−1+Δτk3,Vτj=Vτj−1+Δτ6k1+2k2+2k3+k4,where Gτ,V=A−1MVτ+Lτ+Fτ. The vectors L and F are given by the formulas (24) and (25), respectively, and the matrices A and M are given by the formulas (26) and (27), respectively.


## 3.1. Solving Heat Equation by Compact Finite Difference Method

Combining (1)with formulas (14), (16) and (17), and (20) and (21) generates ℵ×ℵ system of linear ordinary differential equations which can be written in vector matrix form as(23)dVdτ=A−1MVτ+Lτ+Fτ,V0=Φξ,where(i)
V=ϑ1,ϑ2,…,ϑℵT is the approximate solution of heat equation (1)(ii)
Φξ=φξ1,…,φξℵT is the primary dataLet us, respectively, denoteg1τ and g2τ by g1 and g2, and then, the vectors L and F are given by(24)Lτ=1Ds0η2g1+μ1g2h2,m1η2g1+μ1g2,0,…,0,m1μ2g1+η1g2,s0μ2g1+η1g2h2T.And(25)Fτ=fξ1,τ,fξ2,τ,…,fξℵ,τT.A and M are matrices resulting from the sixth-order compact difference scheme defined by(26)A=1l1000…0l1l00…00l1l⋱⋮⋮⋱⋱⋱⋱⋱⋮⋮⋱l1l00…00l1l0…000l11.(27)M=M1,1M1,2M1,3M1,4M1,5M1,6…M1,ℵM2,1M2,2M2,3M2,4M2,5M2,6…M2,ℵm1m2−m3m2m10…00m1m2−m3m2m1⋱⋮⋮⋱⋱⋱⋱⋱⋱00…0m1m2−m3m2m1Mℵ−1,1…Mℵ−1,ℵ−5Mℵ−1,ℵ−4Mℵ−1,ℵ−3Mℵ−1,ℵ−2Mℵ−1,ℵ−1Mℵ−1,ℵMℵ,1…Mℵ,ℵ−5Mℵ,ℵ−4Mℵ,ℵ−3Mℵ,ℵ−2Mℵ,ℵ−1Mℵ,ℵ,where for i=1,3¯(28)M1,2i−1=s2i−1+s04h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M1,2i=s2i+s02h3Dη2f1ξ2i+μ1f2ξ2i.Fori=4,ℵ+1/2−1¯(29)M1,2i−1=s04h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M1,2i=s02h3Dη2f1ξ2i+μ1f2ξ2i,M1,ℵ=s04h3Dη2f1ξℵ+μ1f2ξℵ.The last line for the matrixM is given by: for i=0,2¯(30)Mℵ,ℵ−2i=s2i+1+s04h3Dμ2f1ξℵ−2i+η1f2ξℵ−2i,Mℵ,ℵ−2i−1=s2i+2+s02h3Dμ2f1ξℵ−2i−1+η1f2ξℵ−2i−1,Mℵ,ℵ−6=s04h3Dμ2f1ξℵ−6+η1f2ξℵ−6.Fori=1,ℵ−5/2−1¯,(31)Mℵ,2i−1=s04h3Dμ2f1ξ2i−1+η1f2ξ2i−1,Mℵ,2i=s02h3Dμ2f1ξ2i+η1f2ξ2i.The second line of the matrixM is given by(32)M2,i=m2+m14h3Dη2f1ξi+μ1f2ξi,fori=1,3,M2,2=−m3+m12h3Dη2f1ξ2+μ1f2ξ2,M2,4=m1+m12h3Dη2f1ξ4+μ1f2ξ4.Fori=3,ℵ+1/2−1¯,(33)M2,2i−1=m14h3Dη2f1ξ2i−1+μ1f2ξ2i−1,M2,2i=m12h3Dη2f1ξ2i+μ1f2ξ2i,M2,ℵ=m14h3Dη2f1ξℵ+μ1f2ξℵ.The lineℵ−1 for the matrix M is given by(34)Mℵ−1,i=m2+m14h3Dμ2f1ξi+η1f2ξi,fori=ℵ,ℵ−2,Mℵ−1,ℵ−1=−m3+m12h3Dμ2f1ξℵ−1+η1f2ξℵ−1,Mℵ−1,ℵ−3=m1+m12h3Dμ2f1ξℵ−3+η1f2ξℵ−3,Mℵ−1,ℵ−4=m14h3Dμ2f1ξℵ−4+η1f2ξℵ−4.Fori=1,ℵ−3/2−1¯,(35)Mℵ−1,2i−1=m14h3Dμ2f1ξ2i−1+η1f2ξ2i−1,Mℵ−1,2i=m12h3Dμ2f1ξ2i+η1f2ξ2i.


## 3.2. Resolution of Ordinary Differential System (23)

In order to solve the obtained initial-value problem for a system ofN linear ordinary differential equations (23), we consider the single-step difference schemes for calculating the numerical values of the solution. Let uniform grid space for interval 0,T given by τ0=0 and τj=jΔτ for all j=0,1,…,ℵτ, where Δτ>0 is constant. This difference scheme is a gradual method; we start from the given approximate values Vτj−1 and proceed stepwise, calculating approximate values Vτj at the mesh points. The system (23) is solved using the explicit fourth-order Runge–Kutta scheme which is formulated by(36)k1=Gτj−1,Vτj−1,k2=Gτj−1+12Δτ,Vτj−1+12Δτk1,k3=Gτj−1+12Δτ,Vτj−1+12Δτk2,k4=Gτj−1+Δτ,Vτj−1+Δτk3,Vτj=Vτj−1+Δτ6k1+2k2+2k3+k4,where Gτ,V=A−1MVτ+Lτ+Fτ. The vectors L and F are given by the formulas (24) and (25), respectively, and the matrices A and M are given by the formulas (26) and (27), respectively.


## 4. Numerical Tests

Several numerical tests are presented and discussed in this section in order to illustrate the effectiveness and performance of the proposed method. This section deals with some examples taken from literature to which we compare our approximate solutions.First, in order to compare the exact solution noteduexact and the approximate solution V, we define the absolute error by(37)E=w˜exact−V.The maximal absolute error by(38)Emax=maxi,jw˜i,jexact−Vi,j;i=0,ℵ+1¯andj=0,ℵτ¯.And the relative error by(39)Ere=w˜exact−Vw˜exact.Also, to prove the convergence of the method,O is calculated using the following formula:(40)O=logw˜exact−Vhj−logw˜exact−Vhj+1loghj−loghj+1.The results are obtained for spatial steph=0.1,0.05,0.025,0.0125,0.005 and for temporal step Δτ=0.4h2. The rate of convergence is measured by the value of O, and CPU time is calculated in seconds. All computations and the generation of figures are carried out by MATLAB R2020a on an hp Pavilion, Intel(R) Core(TM) i3-3217U CPU@ 1.80 GHz.Example 1.
We consider equations (1)–(3) with(41)fξ,τ=−2ξ2+τ+1τ+13,φξ=ξ2,f1ξ=f2ξ=ξ,g1τ=−14τ+12,g2τ=34τ+12.
And the exact solution is given by(42)w˜ξ,τ=ξ2τ+12.
In Table1, we display the comparison of the relative errors between our method and the relative errors from [34, 35] of w˜0.6,1 for different values of spatial step h. In Table 2, we exhibit the values of exact and approximate solution for ξ∈0,1 at time τ=1. The last column gives the error committed in each node for h=0.025. The results prove the positive agreement between the exact and approximate solution. Figures 1 and 2 show that numerical solutions reproduce satisfactorily the behaviour of the exact solution of the problem at time τ=1 and h=0.025 and at node ξ=1 and h=0.05, respectively.Table 1
Comparison of relative error forw˜0.6,1 for Example 1.
h
E
r
e in [34]
E
r
e in [35]
PresentEre
0.05
9.95797
e
−
08
2.7
e
−
08
1.13679
e
−
09
0.025
6.22376
e
−
09
1.1
e
−
09
1.83824
e
−
11
0.01
1.59335
e
−
10
7.1
e
−
12
7.20411
e
−
13
0.005
—
4.4
e
−
11
2.71157
e
−
12Table 2
Exact and approximate solution values in different nodes on [0, 1] atτ=1 for h=0.025 for Example 1.
ξ
Exact solution
Approximate solution
E
0.0
0.00000000000
e
+
00
−
9.91029307229
e
−
13
9.91
e
−
13
0.1
2.50000000000
e
−
03
2.50000000523
e
−
03
5.23
e
−
12
0.2
1.00000000000
e
−
02
1.00000000024
e
−
02
2.36
e
−
12
0.3
2.25000000000
e
−
02
2.25000000010
e
−
02
9.95
e
−
13
0.4
4.00000000000
e
−
02
4.00000000002
e
−
02
1.60
e
−
13
0.5
6.25000000000
e
−
02
6.24999999994
e
−
02
5.89
e
−
13
0.6
9.00000000000
e
−
02
8.99999999983
e
−
02
1.65
e
−
12
0.7
1.22500000000
e
−
01
1.22499999996
e
−
01
3.61
e
−
12
0.8
1.60000000000
e
−
01
1.59999999993
e
−
01
7.49
e
−
12
0.9
2.02500000000
e
−
01
2.02499999984
−
01
1.61
e
−
11
1.0
2.50000000000
e
−
01
2.49999999999
e
−
01
9.91
e
−
13Figure 1
Comparison of the exact solution and its approximation in the (a) and absolute error curve in the (b) atτ=1 for h=0.025 for Example 1.
(a)
(b)Figure 2
Comparison of the exact solution and its approximation in the (a) and absolute error curve in the (b) atξ=0.6 for h=0.05 for Example 1.
(a)
(b)Example 2.
We consider equations (1)–(3) with(43)fξ,τ=−exp−τξξ−1+δ26δ2+1+2,φξ=ξξ−1+δ26δ2+1,f1ξ=f2ξ=−δ2,g1τ=g2τ=0.
And the exact solution is given by(44)w˜ξ,τ=exp−τξξ−1+δ26δ2+1;δ=0.12.
To demonstrate the accuracy and efficiency of our method through Example2, we establish in Table 3 the global maximal absolute error when the spatial step takes the values from h=0.1 to h=0.0125. In addition, we give the CPU time in the last column which is compared to results obtained from literature in [17] in Table 4, where we remark that our procedure requires less computational time and achieves more precision than the schemes in [17]. In Table 5, relative errors obtained in [36] are compared with analogous results of our technique for spatial step h=0.05 and h=0.025 for different values of time at node ξ=1.
Based on the results obtained in previous tables, it is clear that our method is more fast and more accurate. The exact and numerical solutions of Example2 are shown in Figure 3 for h=0.05.Table 3
Maximal absolute error and CPU time for Example2.
h, ℵ
E
max
CPU time
0.1,ℵ=9
2.54057
e
−
10
0.10
0.05,ℵ=19
4.33033
e
−
12
0.18
0.025,ℵ=39
7.49123
e
−
14
2.30
0.0125,ℵ=79
6.02296
e
−
15
102.36Table 4
Comparison of the numerical results for Example2.
Methods
E
max
CPU time
Present method forℵ=9
2.54
e
−
10
0.11
Forward Euler forℵ=130 in [17]
1.50
e
−
07
103.75
Backward Euler forℵ=130 in [17]
6.70
e
−
06
81.90
Crank-Nicolson forℵ=130 in [17]
1.30
e
−
06
7.51Table 5
Relative error for different times on [0, 1] atξ=1 for Example 2.
τ
Results in [36]
Present method
Present method
h
=
0.01
ℵ
=
39
h
=
0.025
ℵ
=
19
h
=
0.05
0.1
1.09
e
−
12
5.87
e
−
14
3.67
e
−
12
0.2
1.35
e
−
12
6.87
e
−
14
1.74
e
−
12
0.3
1.35
e
−
12
7.34
e
−
14
4.05
e
−
12
0.4
1.28
e
−
12
7.49
e
−
14
4.10
e
−
12
0.5
1.18
e
−
12
7.59
e
−
14
4.12
e
−
12
0.6
1.07
e
−
12
7.62
e
−
14
4.12
e
−
12
0.7
9.71
e
−
13
7.57
e
−
14
4.12
e
−
12
0.8
8.80
e
−
13
7.65
e
−
14
4.12
e
−
12
0.9
7.96
e
−
13
7.57
e
−
14
4.12
e
−
12
1.0
7.20
e
−
13
7.59
e
−
14
4.12
e
−
12Figure 3
Exact solution in the (a) and its approximation in the (b) forh=0.05 for Example 2.
(a)
(b)Example 3.
We consider equations (1)–(3) with(45)fξ,τ=1+π2expτcosπξφξ=cosπξ,f1ξ=f2ξ=sinπξ,g1τ=expτ,g2τ=−expτ.
And the exact solution is given byw˜ξ,τ=expτcosπξ.
Maximal error, CPU time, and convergence rate for Example3 are tabulated in Table 6 for different values of spatial step h. As expected, the last column indicates that our technique is sixth-order precise, except for the last case h=0.0125. Noting that the average is given by 6.0179 and that for Δτ=0.3h2, we obtained a convergence rate order strictly greater than 6 for all values of ℵ=9.19.39.79. Table 7 exhibits the values of exact and numerical solutions on different nodes in 0,1 at time τ=0.5 and for spatial step h=0.025. It is clear that our method is accurate when compared to the exact solution.
The exact and approximate solutions and absolute error for Example3 are given in Figure 4 when the spatial step h=0.05.Table 6
Maximal absolute error, CPU time, and convergence rate for Example3.
h, ℵ
E
max
CPU time
O
0.1,ℵ=9
7.57818
e
−
06
0.10
—
0.05,ℵ=19
4.49361
e
−
08
0.22
7.3978
0.025,ℵ=39
4.44783
e
−
10
2.34
6.6586
0.0125,ℵ=79
2.78537
e
−
11
103.36
3.9972Table 7
Values of exact and approximate solution for different valuesξ∈0,1 at τ=0.5 and h=0.025 for Example 3.
ξ
Exact solution
Approximate solution
E
0.0
1.64872127070
e
+
00
1.64872127070
e
+
00
1.78
e
−
15
0.1
1.56802710805
e
+
00
1.56802710809
e
+
00
3.36
e
−
11
0.2
1.33384352698
e
+
00
1.33384352705
e
+
00
6.25
e
−
11
0.3
9.69094048058
e
−
01
9.69094048114
e
−
01
5.54
e
−
11
0.4
5.09482891634
e
−
01
5.09482891665
e
−
01
3.14
e
−
11
0.5
1.00955061342
e
−
16
1.25586627697
e
−
15
1.15
e
−
15
0.6
−
5.09482891634
e
−
01
−
5.09482891665
e
−
01
3.14
e
−
11
0.7
−
9.69094048058
e
−
01
−
9.69094048114
e
−
01
5.54
e
−
11
0.8
−
1.33384352698
e
+
00
−
1.33384352705
e
+
00
6.25
e
−
11
0.9
−
1.56802710805
e
+
00
−
1.56802710809
e
+
00
3.36
e
−
11
1.0
−
1.64872127070
e
+
00
−
1.64872127070
e
+
00
1.78
e
−
15Figure 4
Exact solution in the (a) and its approximation in the (b) and absolute error curve in the (c) forh=0.05 for Example 3.
(a)
(b)
(c)Example 4.
We consider equations (1)–(3) with(46)fξ,τ=π2−1sinπξ+cosπξexp−τ,f1ξ=2sinπξ,f2ξ=−2cosπξ,φξ=sinπξ+cosπξ,g1τ=g2τ=0.
And the exact solution is given by(47)w˜ξ,τ=sinπξ+cosπξexp−τ.
Maximal error, CPU time, and convergence rate for Example4 are reported in Table 8 when the spatial length varies from h=0.1 to h=0.0125. It is worth noting from the results obtained in the last column that our technique is seven-order precise for this example. Table 9 depicts the variations of the values of absolute errors at different times in 0,1 and at node ξ=0.25 where we compared our result together with analogous absolute errors in [37], and so we demonstrate that our method is more accurate. The exact and approximate solutions and absolute errors for Example 4 are plotted in Figure 5 when the spatial step h=0.05.Table 8
Maximal absolute error, CPU time, and convergence rate for Example4.
h, ℵ
E
max
CPU time
O
0.1,ℵ=9
8.91994
e
−
06
0.10
;—
0.05,ℵ=19
3.08863
e
−
08
0.22
8.1739
0.025,ℵ=39
1.82998
e
−
10
2.39
7.3990
0.0125,ℵ=79
1.42120
e
−
12
105.73
7.0086Table 9
Absolute error for different times on [0, 1] atξ=0.25 for Example 4.
τ
Results in [37] ℵ=17
Present methodℵ=15
0.1
1.0900
e
−
06
1.35016
e
−
07
0.2
3.0482
e
−
06
1.27234
e
−
07
0.3
9.0128
e
−
06
1.12787
e
−
07
0.4
1.5694
e
−
06
1.00692
e
−
07
0.5
2.0250
e
−
06
9.05450
e
−
08
0.6
2.1913
e
−
06
8.17115
e
−
08
0.7
2.0822
e
−
06
7.38540
e
−
08
0.8
1.6813
e
−
06
6.67954
e
−
08
0.9
1.0599
e
−
06
6.04276
e
−
08
1.0
3.3227
e
−
06
5.46729
e
−
08Figure 5
Exact solution in the (a) and its approximation in the (b) and absolute error curve in the (c) forh=0.05 for Example 4.
(a)
(b)
(c)


## 5. Conclusion

In this paper, the sixth-order compact finite difference scheme is applied to the one-dimensional heat equation with boundary integral conditions, replacing the standard boundary condition. We construct the numerical approach in Section2 and Section 3 in two steps. First, a sixth-order compact finite difference formula is presented in order to obtain a numerical solution to the nonhomogeneous diffusion equation with NCBC. In addition, Simpson’s one-third rule formula is used to meet the difficulty of the integral boundary conditions, while also eliminating additional variables to obtain a square algebraic system. In the second step, the algebraic system obtained is solved using the explicit fourth-order Runge–Kutta scheme, which is unconditionally stable. Overall, the numerical results in Table 2 for Example 1, Table 7 for Example 3, Figure 3 for Example 2, Figure 4 for Example 3, and Figure 5 for Example 4 show that the developed technique gives results much closer to the exact solution, and it demonstrates their high accuracy and efficiency. The strength of this method lies in the fact that it converges fast and requires less computational time, as illustrated in the CPU time calculated in Examples 3 and 4. Moreover, compared to other techniques in the literature, like in Table 4 for Example 2, the proposed method requires less computational time than forward Euler, backward Euler, and Crank-Nicolson methods. Thus, this shows its superiority over other techniques. From the above numerical experiments, we can easily see that the method proposed in this paper obtains Oh6, except for some cases when the spatial step h is very small. The described method can be used to solve many time-dependent problems of different partial differential equations with integral types of boundary conditions in an efficient way.



---

*Source: 8224959-2022-12-30.xml*