Numerical solution for a class of parabolic integro-differential equations subject to integral boundary conditions

Many physical phenomena can be modelled through nonlocal boundary value problems whose boundary conditions involve integral terms. In this work we propose a numerical algorithm, by combining second-order Crank–Nicolson schema for the temporal discretization and Legendre–Chebyshev pseudo-spectral method (LC–PSM) for the space discretization, to solve a class of parabolic integrodifferential equations subject to nonlocal boundary conditions. The approach proposed in this paper is based on Galerkin formulation and Legendre polynomials. Results on stability and convergence are established. Numerical tests are presented to support theoretical results and to demonstrate the accuracy and effectiveness of the proposed method


Introduction
In the last decades, the theory of integrodifferential equations has been extensively investigated by many researchers, and it has become a very active research area. The study of this class of equations ranges from the theoretical aspects of solvability and well-posedness to the analytic and numerical methods for obtaining solutions. A strong motivation for studying integrodifferential equations of PDEs type comes from the fact that they could serve as mathematical models for many problems in physics, mechanics, biology and other fields of sciences.
In this work, we are concerned with the numerical solution of the following parabolic integrodifferential equation: (2.1) where the functional K(·, ·) is defined as follows: Here and in what follows, we use the notation (·, ·) to denote the L 2 -inner product and · for the induced norm on the space L 2 ( ). Denote by H m ( ) the standard Sobolev space with norm and semi-norm denoted by · m and | · | m , respectively. Solvability of the above variational problem is addressed in the following theorem [10].

Space discretization: LC-PSM
Let P N ( ) be the space consisting of all algebraic polynomials of degree at most N and denote by I C N : Based on the above weak formulation, we pose the semi-discrete Legendre-Chebyshev Galerkin schema as: Let L k be the kth degree Legendre polynomial defined by the following three-term recurrence formula: We recall that the set of Legendre polynomials is mutually orthogonal in L 2 ( ), namely Let N be a positive integer, we define [5] The following lemma is the key technique in our algorithm.
Lemma 2.2 [22] For two integer j, k ∈ N, let us denote, and Thanks to linear algebra arguments on can easily prove that Consequently, the numerical solution u N of (2.3) can be expanded in terms of (ϕ k ) N k=0 with time-dependent coefficients, namely (2.5) Inserting (2.5) into (2.3) and taking v = ϕ j , 0 ≤ j ≤ N , we obtain the following system of ODEs Then, the initial value problem (2.6) and (2.7) can be written in matrix formulation as follows: The coefficients m jk and p jk are already determined in Lemma (2.2). For the matrix Q, one can uses the values of φ j (±1) to determinate its entries. In fact, since φ j (±1) = 0 for 0 ≤ j ≤ N − 2, hence Q is almost-null matrix except the two last rows whose entries

Fully-discretization schema
For time advancing, we use the second-order Crank-Nicolson scheme to discretize the differential system (2.8). For a given positive integer M, we define the 3) leads to the following recurrent algebraic system The above algebraic system can be solved easily using either direct or iterative methods. As a choice, on can use QR factorization method, given its accurate results and ease of implementation.

Error analysis
In this section, we derive L 2 -error estimate for the error e N (t) = u N (t) − u(t). For this purpose, we first, in the next subsection, recall a sequence of lemmas that will be needed to perform the error analysis.

Preliminaries
Now, we introduce two projection operators and their approximation properties. First, let P N : L 2 ( ) → P N ( ) be the L 2 -orthogonal projection, namely We also define the operator P 1 N : From the definition of P 1 N , one can obtain Next, we give the approximation property of the projection operator P 1 N and the interpolation operator I C N . Lemma 3.1 [14] If v ∈ H r ( ) with r ≥ 1, then the following estimate holds where C > 0 is a positive constant independent on N .
Lemma 3.2 [14] Let v ∈ H 1 ( ), there exists a positive constant C independent on N such that where C > 0 is a positive constant independent on N . Remark 3.3 Under the same assumptions of Lemma (3.2), we can obtain using approximation property (3.3) the following inequality Now, we derive a basic estimate that will be used later in our proofs.

Error estimates
In this subsection, we consider the stability and convergence of the semi-discrete approximation (2.3). We first state a Gronwall-type inequality that will be used in the proof of our main results.

Lemma 3.5 Let E(t) and H (t) be two non-negative integrable functions on
where C 1 , C 2 ∈ R + . Then there exists C > 0 such that Proof For a non-negative function E(t), we perform a permutation of variables to obtain: Hence, inequality (3.7) of Lemma (3.5) becomes Now, applying the standard Gronwall inequality yields the desired estimate (3.8).
Theorem 3.6 Let u 0 ∈ H 1 ( ) and f ∈ C 1 0, T ; H 1 ( ) . Then the solution u N (t) of (2.3) satisfies (3.10) We have to estimate the terms on the right-hand side of (3.10). For the first term I 1 , we use the hypothesis (1.4) and then apply Cauchy and Young inequalities. (3.11) Next, combining Cauchy and Young inequalities with approximation property (3.5) to estimate I 2 . (3.12) The estimate of I 3 is an immediate consequence of Lemma (3.6), namely Putting things together and choosing 0 < ε < 1 yields (3.14) Integrating both sides of (3.14) form 0 to t, we obtain where Thanks to the Gronwall-type inequality (3.5), we get Lemma 3.7 Assume that u ∈ C 1 (0, T ; H r ( )) , r ≥ 2. Then the following estimate holds where C > 0 is a positive constant independent on N .
Proof From (2.1), (2.3) and (3.1) we know that for a fixed t ∈ J the θ N (t) satisfies for all v ∈ P N ( ) the following error equation: Now, we estimate the terms on the right hand-side of inequality (3.17) using a standard procedure. For the term I 1 , we apply Cauchy and Young inequalities and take into account (1.4), In a similar manner, we can obtain for I 2 By virtue of approximation property (3.2), we bound I 2 as follows For the term I 3 (3.21) Similarly, To estimate of the term I 5 we use Lemma (3.4). Setting w = θ N (t) + ρ N (t) and v = θ N (t) in (3.6) yields using the triangular inequality hence, due to Lemma (3.1), on can obtain, In virtue of above estimates, the inequality (3.18) becomes By taking ε sufficiently small and integrating (3.26) over (0, t), we obtain where Gronwall-type inequality (3.5) implies Take into account, for all 0 < t ≤ T , which is the desired result. Now, we are in position to state our main result concerning the convergence of the semi-discrete approximation (2.3).

Theorem 3.8 Let u(t) and u N (t)
be the solution of (2.1) and (2.3), respectively. If u ∈ C 1 (0, T ; H r ( )) with r ≥ 1, then the following error estimate holds, where C > 0 is a positive constant independent on N .
Proof Using triangular inequality, we have By the aid Lemmas (3.2) and (3.7), for all t ∈ J we obtain This completes the proof.

Numerical experiments
In this section, we carry out several numerical experiments to verify the efficiency and accuracy of the proposed (LC-PSM), and we will compare our results against results obtained using other methods. Example 4.1 In this first test problem, the following parabolic integrodifferential equation is considered The exact solution to the above integrodifferential problem is given as Figure 1 presents the computational results obtained by applying (LC-PSM) to the above test problem, where the profiles of exact and approximate solutions as well as the absolute error are plotted. From the numerical results illustrated in Fig. 1, one can observe that the approximate solution shows a great agreement with the exact solution, which confirms that (LG-PSM) yields a very accurate an efficient numerical method for the numerical resolution of nonlocal boundary value problems of integrodifferential parabolic type.
For comparison purposes, in Tables 1 and 2 we compared our computational results with the results obtained in [10]. Obviously, the proposed (LC-PSM) in this paper gives more accurate solutions with less CPU time than the finite difference schema used in mentioned reference.

Example 4.2
To examine the spatial discretization, we take in this example a test problem that has an analytic solution with limited regularity. Let us consider the following problem: , The exact solution is given as the following: We first choose a step time small enough so that the error of the temporal discretization can be eliminated, and make the polynomial degree N varies. Table 3 shows the error in L 2 and L ∞ -norms at a selected point t = 1 and by going through each line one can observe an increasing accuracy until the error of the temporal discretization becomes dominant.
To examine the theoretical result, we plot in Fig. 2 the decay rates of error in L 2 -norm versus N in a logscale and the lines of decay rates N −2 and N −4 . As expected, L 2 -error of (LC-PSM) for the solved problem in this example has a rate of convergence between N −3 and N −4 , which supports the results established in Theorem (3.8) since u ∈ H 3 ( ) and u / ∈ H 4 ( )

Conclusions
In this paper, we are concerned in the implement and analysis of the spectral method to solve a class of integrodifferential parabolic equations subject to nonlocal boundary conditions of Neumann-type. We combined the Legendre spectral method based on Galerkin formulation to discretize the problem in the spatial direction and the second-order Crank-Nicolson finite difference schema for the temporal discretization. Rigorous error analysis has been carried out in L 2 -norm for the proposed method, and the computational results of numerical examples have supported the theoretical results. Moreover, a comparison with fully finite-difference schema clearly shows that the presented method is computationally superior with less required CPU time. It should be noted that other high-order methods can be used for time integration to improve the accuracy of the fully discretization. Convergence and stability of such combinations are still undiscussed.
In future works, we plan to investigate how to implement space-time spectral method for the resolution of this class and other challenging models, such as nonlocal boundary value problems in the two-dimensional case and fractional integrodifferential problems.