Error Analysis of Legendre-Galerkin Spectral Method for a Parabolic Equation with Dirichlet-Type Non-Local Boundary Conditions

. An eﬃcient Legendre-Galerkin spectral method and its error analysis for a one-dimensional parabolic equation with Dirichlet-type non-local boundary conditions are presented in this paper. The spatial discretization is based on Galerkin formulation and the Legendre orthogonal polynomials, while the time derivative is discretized by using the symmetric Euler ﬁnite diﬀerence schema. The stability and convergence of the semi-discrete spectral approximation are rigorously set up by following a novel approach to overcome diﬃculties caused by the non-locality of the boundary condition. Several numerical tests are included to conﬁrm the eﬃcacy of the proposed method and to support the theoretical results.


Introduction
Initial boundary value problems for partial differential equations with non-local boundary conditions NLBCs have been extensively investigated in many papers and textbooks. The main reason of such enormous attention for this kind of problems is that some problems in physics, chemistry and many other fields of the sciences can be effectively approached using models based on non-classical boundary value problems for PDEs [6,12,13,15,19,25]. This paper is devoted to develop and analyse the implementation of Legendre-Galerkin spectral method to non-local boundary value problem for a linear parabolic equation. Our problem is given as follows ∂u ∂t (x, t) − ∂ 2 u ∂x 2 (x, t) = f (x, t), (x, t) ∈ Λ × (0, T ], u(x, 0) = u 0 (x), x ∈Λ u(x, t)K 2 (x)dx, (1.2) where Λ = (−1, 1) and T stands for a final time. The integral kernels K 1 and K 2 , the source term f and the initial data u 0 are given smooth functions.
Recently, researchers have shown an increased interest in analysing and developing computational methods for the numerical solution of nonlocal boundary value problems for partial differential equations, including the diffusion equation subject to nonlocal boundary condition of the form (1.2). Considerable number of existing papers that dealt with the numerical solution for this kind of problems are mainly based on the finite difference schemas, in which the stability, convergence and other numerical aspects were considered (see e.g., [7,8,9,21,23] and the references therein). By way of example, Jakubėlienė et al. [16] have studied the stability and convergence of difference schemes for approximating a class of semilinear parabolic equations with one nonlocal condition. In the same vein, discussions on the impact of complex coefficients on the stability of difference schemes were reported in [22]. In a recent study investigating a semi-implicit difference scheme for a two-dimensional parabolic equation with NLBCs, the authors [10] applied a methodology based on the properties of M-matrices to examine the stability of the proposed finite difference scheme. However, solving partial differential equations with either local or nonlocal boundary conditions by mesh-dependent schemas such as finite difference methods and finite elements methods is computationally very expensive, as discretization for both space and time variables is required.
In the same context, the familiar operational matrices approaches have been used by many authors to handle with several problems that fall in this category. The authors of [30] adopted a new technique based on operational matrices of Bernstein polynomials to approximate the one-dimensional diffusion equation with an integral condition. Recently, Borhanifar et al. [4] presented an effective algorithm based on operational matrix formulation to solve nonlinear reactiondiffusion equations with mixed nonlocal boundary conditions, despite its excellent computational results for smooth problems, it seems inappropriate to some special kind of PDEs such as oscillatory problems. Other various methods have been also applied to the nonlocal boundary value problem (1.1)-(1.2) and similar problems, such as radial basis functions method [17,29], reproducing kernel procedure [20,32], and so on.
Compared to the aforementioned numerical methods, Galerkin spectral methods [5,27], whose accuracy and robustness were widely recognized, can provide excellent computational results that enjoy superior accuracy with only moderate discretization resources. Despises this fact, there is only a limited number of papers that touched upon the implementation, in particular, the analysis of spectral methods for PDEs with non-local boundary conditions, which is basically due to the lack of a theoretical framework corresponding to the non-locality of the boundary conditions.
To our best knowledge, there is no research available in the literature devoted to analysis the implementation of spectral methods to nonlocal boundary value problems for parabolic equations. Motivated by this fact, we mainly aim in this paper to provide a suitable approach to solve the one-dimensional parabolic equation (1.1) subject to non-local boundary conditions (1.2) by a spectral method with efficient implementation and exponential rate of convergence as in the spectral methods for problems with classical boundary conditions. We emphasize the fact that the parabolic equation with boundary conditions (1.2) considered in this work is not only interesting in its theoretical and practical applications, but also in the methodology proposed to handle with this problem, which constitutes a flexible approach that can be extended for solving more general PDEs subject to NLBCs (see e.g., [2,11,26]).
The objective of this paper is twofold. First, to construct an appropriate spectral method for the parabolic equation (1.1) subject to boundary conditions (1.2). Second, to carry out some error analysis for the proposed method. A drawback of using spectral Galerkin formulation in a standard way for nonlocal boundary value problems is, of course, the lack of an appropriate basis whose elements satisfying the integral conditions at the boundary. The key technique to overcome this obstacle is to use a non-classical weak formulation that treats the main equation and the boundary conditions separately. As we will see in the paper, the use of such formulation also allows us to perform the error analysis.
The rest of the paper is organized as the following. First, in the next section, we briefly introduce the Legendre-Galerkin spectral method for solving the problem (1.1)-(1.2) and describe the way to implement the proposed method. Section 3 is devoted to recollect some fundamental notions and results of spectral approximation needed for the error analysis. In Section 4, we provide error bounds in L 2 and H 1 -norms for the semi-discrete Legendre-Galerkin spectral method. Then, in Section 5, we examine the accuracy and robustness of our method by extensive numerical tests. We conclude the paper with some remarks on the main features of the method presented in previous sections and highlight some possible extensions of our method.

Legendre-Galerkin spectral method
We aim in this section to give a brief description of the implementation of the Legendre-Galerkin spectral method (LG-SM) to solve the problem under consideration. Through this paper, (·, ·) and · denote scalar product and its associated norm in L 2 (Λ). We also use H m (Λ)(m ≥ 1) to denote the standard Sobolev space endowed with the norm . m and the semi-norm | · | m , where We define the space Let P N (Λ) stands for the space of all algebraic polynomials of degree at most N , we also define the space P 0 N (Λ) = P N (Λ) ∩ H 1 0 (Λ). In order to derive the form of our approximation schema, we first reformulate problem (1.1)-(1.2) into its weak formulation, namely, find u(t) ∈ H 1 (Λ) such that For the solvability of the above weak formulation, Slodicka [28] proved, by the use of Rothe method, the existence and uniqueness of the solution.
The conventional Legendre-Galerkin semi-discrete approximation reads as: where i C N is the operator of interpolation at the Chebyshev-Gauss-Lobatto nodes ξ j = cos ( jπ N ), 0 ≤ j ≤ N . Let us denote L k (x) the k-th degree Legendre polynomial. The set of Legendre polynomials {L k } ∞ k=0 forms an orthogonal basis for the space L 2 (Λ), namely, Otherwise, any function u in L 2 (Λ) can be expensed in terms of L k , Let N be a positive integer, we define Obviously, the set {φ k } N k=0 consists of N + 1 linearly independent elements, therefore form a basis function for P N (Λ). We set Upon substituting (2.3) in the boundary conditions, we get two supplementary equations, namely, For the time advancing, we use the second-order symmetric Euler method for system (2.4) with boundary conditions (2.5). Let M be positive integer and define a step time ∆t = T M . Let t i = i∆t(0 ≤ i ≤ M ), we denote byû i k the approximation ofû N (t k ). At each time level, we need to solve the following algebraic linear system, By the aid of the properties of Legendre polynomials, on can easily determinate the values of the quantities m j,k and p j,k , as follows [14,31]

Error analysis
This section is devoted to set up the convergence and stability results for semidiscrete Legendre-Galerkin spectral method, we also aim to show that the proposed method enjoys the spectral accuracy. In what follows, C will denote a generic positive constant not depending on N and it may take different values each time we use it.

Preliminaries
In order to carry out error estimates for the schema (2.2), we introduce three projection operators with their corresponding approximation properties. We define the L 2 -orthogonal projection P N : L 2 (Λ) → P N (Λ) as the following be the orthogonal projection operator defined as the following It is readily verified that Next, we give approximations proprieties of P N and P 1,0 N . Lemma 1 [ [3], [5]]. Let r and s be two non-negative real numbers with r ≤ s. For any u ∈ H s (Λ), the following estimates hold where C > 0 is positive constant depending only on s. [5]]. Let r and s be two non-negative real numbers, with the assumption where C > 0 is positive constant depending only on s.
In the context of the numerical resolution of non-local boundary value problems, we generally deal with the approximation of functions which do not vanish at the boundary points. For this purpose, we need to define an appropriate projection operator which preserve the function values at the boundary.
Following the same idea as in [5, page 288], for any function u in H 1 (Λ), we first introduce the polynomial and then we associate for u the functionū := u−K[u]. Since K[u](−1) = u(−1) and K[u](1) = u(1), thusū ∈ H 1 0 (Λ). This allows us to introduce the following operator P 1,b N : H 1 (Λ) → P N given by We give the approximation properties for the operator P 1,b N by stating the following lemma Lemma 3. For any function u ∈ H 1 (Λ), . Moreover, for any real numbers r and s with 0 ≤ r ≤ 1 ≤ s. If u ∈ H s (Λ) then the following estimate holds

4)
where C > 0 is a positive constant depending only on s.
We conclude this section by giving some approximation properties for the operator of interpolation i C N . Lemma 4. For any u ∈ H 1 (Λ) the following estimate holds where C > 0 is a positive constant independent on N .

Stability and convergence
This section is devoted to analysing the stability and the convergence rate of the semi-discrete scheme (2.2). We first give a theoretical analysis for the stability by stating two lemmas hereafter. Before this, we derive some basic inequalities which are needed in our proofs in the sequel. For convenience, we denote p(x) = 1 2 (1 − x). Using Cauchy and triangle inequalities, (3.7) Analogously, we can bound the term ∂ x K[u N (t)] as the following where: α K = √ 6 Math. Model. Anal., 26 (2):287-303, 2021.
Lemma 5. Let u 0 ∈ H 1 (Λ), ∂ t u(t) ∈ C 1 (0, T ; H 1 (Λ)) and f ∈ C 1 (0, T ; L 2 (Λ)). If α K < 1, then there exists a positive constant C independent on N such that the solution u N (t) of the semi-discrete approximation (2.2) satisfies the following estimate where We shall bound the terms I 1 , I 2 and I 3 in a standard way through both Cauchy and ε-Young inequalities using basic inequalities (3.7) and (3.8).
We begin with I 1 .
Next, we estimate I 2 The term I 3 can be controlled in the same manner. (3.9) Putting things together and fixing 0 < ε < 1 − α K , we get Integrating inequality (3.10) over (0, T ) and taking into account approximation property, we obtain By the aid of Gronwall inequality, we can get the desired result.
Lemma 6. If u 0 ∈ H 1 (Λ) and f ∈ C 1 (0, T ; L 2 (Λ)), then there exists a positive constant C independent on N , such that the solution u N (t) of the semi-discrete approximation (2.2) satisfies the following estimate where Estimates for the terms on right hand-side of (3.9) can be derived essentially in the same manner as in the proof of the previous lemma. As a matter of fact, using Cauchy and Young inequalities together with basic inequalities (3.7) and (3.8), we can get Putting things together and fixing ε sufficiently small, we get The integration of (3.11) combined with (5) yields We see that u N (0) = i C N u 0 = i C N u 0 − u 0 + u 0 and by the use of (3.6), it follows that u N (0) ≤ C u 0 1 . Therefore, the proof is completed by applying Gronwall inequality. Now we turn to the convergence analysis. We derive error estimate in L 2 and H 1 -norms of the error E N (t) between the semi-discrete solution u N (t) of the schema (2.2) and the exact solution of the problem (2.1). To this end, for any According to approximations properties (1) and (2), we only need to bound θ N (t). Before stating our main results, we first need to derive some auxiliary estimates. By a technical reduction, one can get Using Cauchy and triangle inequalities, we derive from the above identity the following basic estimates , , (3.12) where α K and β K stand for the same constants in (3.7) and (3.8).
Lemma 7. Let r ≥ 1 be a positive real number. Assume that u ∈ H 1 (0, T ; H r (Λ)) and f ∈ H 1 (0, T ; H r (Λ)). If α K < 1, then the following estimate holds Proof. For an arbitrary function ϕ ∈ P 0 N (Λ), we difference (2.1) and (2.2) and where Estimates for the terms on the right hand-side of (3.13) are quite similar to those which appear in the proofs of previous lemmas. In analogue manner, by the aid of Cauchy and ε-Young inequalities besides basic estimates (3.12) we can get Combining the above estimates with (3.13), to get By choosing 0 < ε < 1 − α K and utilizing Gronwall inequality, we get (3.14) Meanwhile, we can see that r . Therefore, inequality (3.14) takes the form which is the desired result.
Lemma 8. Let r ≥ 1 be a positive real number. Assume that u ∈ H 1 (0, T ; H r (Λ)) and f ∈ C 1 (0, T ; H r (Λ)) for certain r ≥ 1. Then the following estimate holds Proof. The proof of Lemma (8) can be done exactly in the same manner as the one of the previous Lemma.
Next, on the basis of Lemmas (7) and (8) we state the following result directly without proof Lemma 9. Let r ≥ 1 be a positive real number. Assume that u ∈ H 1 (0, T ; H r (Λ)) and f ∈ C 1 (0, T ; H r (Λ)). Then the following estimate holds θ N (t) 2 1 ≤ CN −2r , where C is a positive constant not depended on N . Now we are in position to provide error bounds for the semi-discrete approximation (2.2) by stating the following theorem. Theorem 1. Let r ≥ 1 be a positive real number. Assume that u ∈ H 1 (0, T ; H r (Λ)), f ∈ C 1 (0, T ; H r (Λ)) and α K < 1. Then the following estimates hold

16)
where C is a positive constant not depended on N .
Proof. By definition, Applying triangular inequality and using Lemmas (8) and (3), we get the L 2error estimate, Following an analogue manner we can prove the H 1 -error estimate by the aid of Lemma (9).

Numerical experiment
In this section, we present several numerical tests using Legendre-Galerkin spectral method described in previous sections. During the process, we computed the integral terms using the Clenshaw-Curtis quadrature for high accuracy. Meanwhile, the linear algebraic systems resulted from the full-discrete approximation had solved by applying a direct method based on LU decomposition.
In this first example, we begin our computational investigation by testing the accuracy to demonstrate the effectiveness of our method. Set N = 16 and ∆t = 10 −3 . Figure 1 displays the profiles of exact and approximate solutions, and also shows the absolute error function E N = |u − u N |. The results show an excellent agreement between the exact and approximate solutions, which confirm the accuracy of (LG-SM) method.  We turn to the convergence behaviour of the approximate solution. For the spatial discretization, we fix the step time ∆t sufficiently small, say ∆t = 10 −5 , so that the error of temporal discretization is negligible, and vary N ≤ 20.
The error u − u N in L ∞ and L 2 -norms at two selected points t = 1 and t = T are reported in Tables 1 and 2. According to the computational results, the error shows an exponential decay, which demonstrates the spectral accuracy property. Next, we check the temporal discretization. By choosing N big enough, so that, the spatial discretization error is negligible, and let ∆t varies. The results are presented in Table 3. As expected, the temporal accuracy is of order 2, which is consistent with theoretical results.

Example 2.
As for all other spectral methods, the convergence rate of the present spectral method is depending on the regularity of the solution. In this example, we consider a problem with a non-smooth solution to examine the validity of the error estimates derived above. Let us consider the following exact solution with limited regularity of the Equation (1.1) subject to the following boundary conditions In Figure 2, we plot error decay rates in L 2 and L ∞ -norms versus the polynomial degree in log-scale with a fixed step time ∆t = 10 −4 . The lines of decay rates N −3 and N −5 are also plotted for a close comparison. It is obvious that the L 2 and L ∞ -errors decay with a rate between N −3 and N −5 , which is match with error estimate (3.16). This computational result seems better than the theoretical one, as the solution belongs at most to H 2 (Λ).
where the kernels in the nonlocal boundary conditions (1.2) are chosen to be The exact solution of this problem is given as, u(x, t) = −e t x(x − 1) + δ 6(1 + δ) .
For this problem, we applied Legendre-Galerkin spectral method described in the previous sections by taking N = 4 and ∆t = 10 −5 and compare the computational results for δ = 0.0144 with the results obtained in [1,4,24] (see Tables  4 and 5). Obviously, the computational results of LG-SM are better than the obtained ones using the algorithms presented in [1,4,24] even for small polynomial degree N , this may relax the storage limit and decrease the computational cost. Moreover, we can obtain more accurate results if a high-order schema is employed for the time-discretization. Table 5. Absolute errors of Example 3 using methods in [4,24] and LG-SM for t = 0.5. x Method in [24] Method in [4] LG-SM N = 10 N = 8 N = 4

Conclusions
In the present paper, Legendre-Galerkin spectral method for the diffusion equation with non-local boundary conditions is proposed. Based on the weak formulation of the problem under consideration, we were able to implement and analyse the proposed spectral method. Our approach enjoys the spectral accuracy property, which has proven by deriving suitable error estimates, and illustrated by several numerical examples. This work is a primary attempt in developing spectral methods for parabolic equations with non-local boundary conditions, in which, we only considered a one-dimensional linear parabolic equation. The approach developed in this work could serve as a base for the analysis and implementation of spectral methods for more general PDEs with analogous boundary conditions, such as semilinear and nonlinear parabolic equations. In further works, we plan to extend the present approach to investigate space-time spectral methods for solving such classes of PDEs.