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

#### 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:where is a sufficiently differentiable function in time and space, , and subject to the initial conditionand the NCBCwhere , , , , 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 , with boundary condition and subject to the energy specification 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., is replaced by with homogeneous initial condition and the boundary condition as well as the NCBC(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 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., . 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 with nonlocal boundary condition (5), such that vary from to , the initial condition (2), and the boundary condition(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 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 Section 2 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 interval into subintervals having length , and the time interval is divided into subintervals having length . The grid points are detailed by

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

For , the scheme is sixth-order accurate, and in this case, , , and the truncation error is .

By simplifying the relation (13), we obtainwhere

In order to get the same accuracy and dominant error term, special formulas are needed for and . When , we use

When , we usewhere the coefficients can be found by matching the Taylor series expansions up to the order of , which gives us a linear system in which the solution is given by

All formulas are . More details are available in [24–26]. When unknowns and are managed by the boundary conditions of the considered problem, we obtain a system of equations with unknowns , , .

#### 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 constructs a system of linear equations with unknowns , , , . Simpson’s one-third rule formula is used to approximate the integrals in (3) and thus eliminate and , to obtain a system of ordinary differential linear equations in unknowns , , .

So, we havewhere

##### 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 aswhere(i) is the approximate solution of heat equation (1)(ii) is the primary data

Let us, respectively, denote and by and , and then, the vectors and are given by

And

and are matrices resulting from the sixth-order compact difference scheme defined bywhere for

For

The last line for the matrix is given by: for

For ,

The second line of the matrix is given by

For ,

The line for the matrix is given by

For ,

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

In order to solve the obtained initial-value problem for a system of 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 given by and for all , where is constant. This difference scheme is a gradual method; we start from the given approximate values and proceed stepwise, calculating approximate values at the mesh points. The system (23) is solved using the explicit fourth-order Runge–Kutta scheme which is formulated bywhere . The vectors and are given by the formulas (24) and (25), respectively, and the matrices and 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 noted and the approximate solution , we define the absolute error by

The maximal absolute error by

And the relative error by

Also, to prove the convergence of the method, is calculated using the following formula:

The results are obtained for spatial step and for temporal step . The rate of convergence is measured by the value of , 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) withAnd the exact solution is given byIn Table 1, we display the comparison of the relative errors between our method and the relative errors from [34, 35] of for different values of spatial step . In Table 2, we exhibit the values of exact and approximate solution for at time . The last column gives the error committed in each node for . 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 and and at node and respectively.

**(a)**

**(b)**

**(a)**

**(b)**

*Example 2. *We consider equations (1)–(3) withAnd the exact solution is given byTo demonstrate the accuracy and efficiency of our method through Example 2, we establish in Table 3 the global maximal absolute error when the spatial step takes the values from to . 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 and for different values of time at node .

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 Example 2 are shown in Figure 3 for .

**(a)**

**(b)**

*Example 3. *We consider equations (1)–(3) withAnd the exact solution is given by

Maximal error, CPU time, and convergence rate for Example 3 are tabulated in Table 6 for different values of spatial step . As expected, the last column indicates that our technique is sixth-order precise, except for the last case . Noting that the average is given by and that for , we obtained a convergence rate order strictly greater than 6 for all values of . Table 7 exhibits the values of exact and numerical solutions on different nodes in at time and for spatial step . It is clear that our method is accurate when compared to the exact solution.

The exact and approximate solutions and absolute error for Example 3 are given in Figure 4 when the spatial step .

**(a)**

**(b)**

**(c)**

*Example 4. *We consider equations (1)–(3) withAnd the exact solution is given byMaximal error, CPU time, and convergence rate for Example 4 are reported in Table 8 when the spatial length varies from to . 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 and at node 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 .

**(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 Section 2 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 , except for some cases when the spatial step 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.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that they have no potential conflicts of interest.

#### Acknowledgments

N. Arar and L. Ait Kaki acknowledge support from the Directorate General for Scientific Research and Technological Development (DGRSDT), Algeria.