Rothe–Legendre pseudospectral method for a semilinear pseudoparabolic equation with nonclassical boundary condition

. A semilinear pseudoparabolic equation with nonlocal integral boundary conditions is studied in the present paper. Using Rothe method, which is based on backward Euler ﬁnite-difference schema, we designed a suitable semidiscretization in time to approximate the original problem by a sequence of standard elliptic problems. The questions of convergence of the approximation scheme as well as the existence and uniqueness of the solution are investigated. Moreover, the Legendre pseudospectral method is employed to discretize the time-discrete approximation scheme in the space direction. The main advantage of the proposed approach lies in the fact that the full-discretization schema leads to a symmetric linear algebraic system, which may be useful for theoretical and practical reasons. Finally, numerical experiments are included to illustrate the effectiveness and robustness of the presented algorithm.


Introduction
Numerous scientific and engineering problems can be modelled employing nonlocal partial differential equations. This, in turn, results in a need for examination and investigation of the behavioural and qualitative properties of these models.
In this paper, we concerned by the solvability of a semilinear pseudoparabolic equation accompanied by nonclassical boundary conditions. Let Ω be a bounded subdomain of R d , (d 1) with a Lipschitz continues boundary ∂Ω.
Here I T stands for the time interval [0, T ] with T > 0, and F supposed to be a Lipschitz continuous function. Historically, the first appearance of the term pseudoparabolic equation was in the works of Ting and Showalter [21,26], since then, the pseudoparabolic equation is referred to a partial differential equation with a first-order derivate in time appearing in the highest-order term, in particular, such kind of problems forms a subclass of Sobolev-type equations [25].
This kind of PDEs appears, while studying variety of physical and engineering models, for example, transport phenomena in porous media [2,8,10], the motion of homogeneous fluids in fissured rock [1], moisture transfer in soils [19], the diffusion of imprisoned radiation through a gas [5,11,12], and so on.
Due to its applications in several fields of applied sciences, initial value problems for pseudoparabolic equations subject to classical, as well as nonclassical, boundary conditions are of important practical and theoretical interest and have been investigated in theoretical and numerical frames by several authors. Some of those authors have studied the solvability of a variety of initial boundary value problems for the third-order equations of pseudoparabolic type. For example, Popov [20] investigated the one-dimensional linear pseudoparabolic equations with integral conditions for which he used an approach based on the method of continuation with respect to a parameter and some a priori estimates. The authors in [7] have established, using some appropriate elliptic estimates, the wellposedness of a nonlinear pseudoparabolic equation with two types of nonlocal boundary conditions. Bouziani and Merazga [3] have used backward Euler schema time discretization to prove the existence of the solution to a similar nonlocal problem in nonclassical space. In the same line, Guezane-Lakoud and Belakroum [9] proved the existence and uniqueness of solution in weak sense to an integrodifferential Sobolev-type equation with integral conditions. Recently, Maqbul and Raheem [17] studied by means of Rothe method the questions of the existence and uniqueness of solution for a semilinear pseudoparabolic equation with purely integral conditions. Along a different line, various numerical methods are aroused to solve nonlocal boundary-value problems. Up to now, the work done on the numerical analysis for nonlocal pseudoparabolic equations is mainly based on the finite-difference methods [6,[13][14][15].
The primary objective of this work is to perform a computational study on problem (1)- (3). Applying a similar methodology as in [23,24], we investigate the existence and uniqueness as well as the approximate construction of solutions and its convergence. In short, we will use the semidiscretization in time to approximate the nonlocal problem by a sequence of elliptic boundary-value problems with Dirichlet boundary condition and apply the Legendre pseudospectral method for the spatial discretization. By deriving some a priori estimates we will able to prove the convergence of the approximations sequence to a unique solution.
It must be emphasized that most researches concerning the nonlocal third-order pseudoparabolic equations have mainly considered purely integral conditions. Meanwhile, the nonlocal boundary-value problem that is considered in this work is accompanied by the nonlocal Dirichlet-type boundary condition (2). Moreover, to the authors best knowledge, there are no works deal with the implement and analysis of spectral and pseudospectral methods for the problem under consideration, even not in the case of one-dimensional setting, which is mainly due to nature of the boundary conditions. Therefore, this work is aimed at extending and improving some existing results on the solvability of nonlocal boundary-values problems of pseudoparabolic type as well as presenting an efficient numerical schema based on coupling the Euler finite-difference method for time discretization with Legendre spectral method for space discretization to obtain approximate solutions for this kind of problems.
Besides the introduction, this paper includes three sections. First, we turn the problem under consideration into an appropriate weak formulation, and then backward finitedifference schema is employed to approximate this last problem by standard elliptic problems in Section 2. In the same section, some stability results are established, and finally, prove the convergence of the semidiscrete temporal approximations to the unique solution of the original problem. In Section 3, we briefly described a proposed way to implement Legendre pseudospectral for the discretization in the space direction. Numerical tests are presented in the last section to illustrate the effectiveness of the presented numerical algorithm.

Time discretization
The study in this paper will be done in variational framework. For this purpose, we need to introduce some function spaces and related useful results.
Let us denote by L 2 (Ω) the space of square-integrable functions on Ω endowed with the usual L 2 -scalar product (·, ·) and its corresponding norm · . Also, we denote by H m (Ω) the standard Sobolev space defined by which is a Hilbert space for the usual scalar product where the derivative ∂ α is taken in weak sense, the norm of the space H m (Ω) is denoted by · m . Now, we state a practical result with a crucial role in our proofs later.
Next, we will give the definition of a weak solution of problem (1)-(3) in the following sense.
Definition 1. We say that the function u is a weak solution to the nonlocal problem (1)- (C3) the following identity holds for all φ ∈ H 1 0 (Ω) and a.e. t ∈ I T : In what comes next, we will denote by C, C ε , and ε generic constants, which are independent of any discretization parameters, where C ε = C(1/ε), and ε is arbitrarily small. These constants may change its values from line to another but the meaning will be clear from the context.

Construction of approximate solution and convergence
In order to prove the existence of a weak solution in the sense of Definition (1), we will use Rothe method. The idea behind this approach is to approximate the time-dependent problem by a sequence of elliptic problems that have to be solved successively at each time level. To this end, we first discretize the continuous problem (6) in time. We divide the interval I T into n ∈ N * equidistant subintervals [t i−1 , t i ], i = 1, . . . , n, where t i = iτ and τ = T /n. Moreover, we adopt the following notations [18]: We approximate the variational problem (6) by the sequence of standard elliptic boundaryvalues problems. Using the above notations, the recurrent semidiscrete approximation scheme for i = 1, . . . , n reads as and the boundary condition (2) can be treated in the following way: The well-posedness of problems (7) i -(8) i is addressed in the following lemma Lemma 2. Let F be Lipschitz continuous function, 0 < α, β C, and assume that u 0 ∈ H 2 (Ω). Then, for any i = 1, . . . , n, there exists a unique solution u i to the variational problem (7) Using estimates (4), (5) and Cauchy inequality, we can easily check that A(·, ·) is a coercive continuous bilinear form on H 1 0 (Ω) × H 1 0 (Ω) and F(·, ·) is a linear bounded functional on H 1 0 (Ω). Thus, by applying Lax-Milgram lemma we obtain the existence of a unique solution v i for (9), that ensures the well-posedness of recurrent problem (7) Now, we derive some a priori estimates by stating the following lemma.
Lemma 3. Let the assumptions of Lemma 2 be satisfied. Moreover, assume that α 2 +β 2 < 1. Then there exists a positive constant C > 0 such that Proof. Testing the identity (7) with φ = (δu i − δK i )τ , which belongs to H 1 0 (Ω), and summing up for i = 1, . . . , j, we get We recall the following identity: where (a i ) is a sequence of real numbers. We use this formula for the left-hand side of (10) to obtain L.H.S. of (10) = 1 2 Now we estimate each terms on the right-hand side of equation (10). Following a straightforward procedure by applying Cauchy and ε-Young inequalities together with basic estimates (4), (5), we can obtain Summarizing the above estimates, R.H.S. of (10) Combining (11) and (12) with the choice of suitable ε < 1/2 − (α 2 + β 2 )/2 and using Gronwall lemma, we get the desired result. We also define the piecewise constant in time functionū n : [0, T ] → H 1 (Ω), u n (0) = u 0 ,ū n (t) = u i ∀t ∈ (t i−1 , t i ], 1 i n. Thus, the weak problem (7)- (8) can be written in terms of the above notation as Now, we are about to prove the convergence of the approximate solution u n to a weak solution to the variational problem (6). First, we reformulate the estimate from Lemma 3 as follows: The following lemma is a direct consequence of (14). Proof. It suffices to see that for all t ∈ I T , In view of Lemma (2), we obtain consequently, u n (t) −ū n (t) 1 → 0 as n → ∞, which is the desired result.
Regarding estimate (14), we see that all conditions of [16, Lemma 1.3.13] are met, which allows us to state directly without proof.
Lemma 5. Under the assumptions of Lemma 3, there exists a function u ∈ C(0, T ; H 1 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)) obeying ∂ t u ∈ L 2 (0, T ; H 1 (Ω)) and a subsequence (u n k ) of (u n ) such that: Now, we are ready to state our main result in this section. Theorem 1. The limit function u from Lemma 5 is the unique weak solution to problem.
Proof. Integrating (13) over (0, T ) and using the fact that u n (0) = u 0 , we obtain Make use of Taking the limit for n → ∞ in (15) and using the results of Lemma 3, we obtain Differentiating (16) with respect to the time variable implies that u is a weak solution to variational problem (6). It remains to check the behavior of (u n ) on the boundary ∂Ω. Following a similar argument as in [22], we start from the inequality Using the trace inequality to get Now, let t ∈ [t i−1 , t i ], we have by definitionū n (t)| Γ = Kū (t−τ ) , which yields after performing some arrangement Next, we estimate each term on the right-hand side of the above identity separately: In similar manner, one can obtain Putting things together, Due to relations (17) and (18), we deduce that Then, passing to the limit n → ∞, we obtain by Lemmas 4 and 5 This proves that the limit u is a solution to (6).
Proof. We suppose, for the sake of contradiction, that there are two solution u and v satisfying (6), and we put w = u − v. Subtracting the corresponding weak formulation for u form the one for v and taking the test function φ = ∂ t (w − K w ), we get for all t ∈ I T , Using Cauchy and ε-Young inequalities and Lemma 1, we can bound the right-hand side of the above identity as follows: Using this estimate an fixing ε < 1 − (α + β), identity (19) becomes Integrating the previous inequality over the interval I T and applying Gronwall lemma gives This yields that u(x, t) − v(x, t) = w(x, t) = 0 almost everywhere in Ω × (0, T ), that achieves the desired result.

Full discretization
In this section, we deal with the spatial discretization of (7), which completes the full discretization of the nonlocal problem (1)- (3). To this end, we apply the Legendre pseudospectral method. For the numerical resolution, we keep things simple by considering only the onedimensional setting. For this purpose, the space domain Ω is chosen to be the finite interval (−1, 1). Moreover, the boundary conditions become Now, it is convenient to reduce the recurrent problem with inhomogeneous conditions to an equivalent problem with homogeneous Dirichlet condition. Let us set Let us denote by L k (x) the kth degree Legendre polynomial and P 0 M the set of all algebraic polynomials of degree at most M defined on (−1, 1), which vanish at the boundary. Then the standard Legendre-Chebyshev approximation to (7) is to find for all i = 1, . . . , n, a functionû M i ∈ P 0 M such that starting withû M 0 = I C M u 0 , and the additional functional G(·) collect the terms resulting from the reformulation. The following lemma is the key to the efficiency of our algorithms.
Lemma 6. (See [4,27].) Let us denote Then According to linear algebra arguments, it is clear that the set {φ k } M −2 k=0 spans the space P 0 M , that allows to writeû M i as a linear combination of φ k , namely, Inserting (22) into identity (21) and taking v = φ j , j = 0, . . . , M − 2, leads to the following algebraic linear system: where The linear system (23) can be easily solved since the matrices M and P are banded. We can use a direct or either iterative method to solve it.

Numerical experiments
The aim of this part of paper is to present some numerical experiments that we carried out in order to confirm the reliability and viability of the proposed procedure. As mentioned above, we consider problem (1) with one-dimensional space variable subject to nonlocal boundary (20).
Example 1. For the first example, we seek for the numerical solution to the problem with the exact solution: u * (x, t) = (sin (πx) + cos (πx))e −t . To solve the above problem, we applied the Rothe-Legendre pseudospectral method described in Sections 2 and 3. In Fig. 1, we plot the absolute error as a function of x and t with fixed polynomial degree N = 12 and various step time τ . The results shows that the approximate and exact solutions are in good agreement, that demonstrates the effectiveness of the proposed algorithm.
To investigate the convergence rates of temporal discretization, we choose the polynomial degree M big enough, say M = 16, so that the spatial discretization error is negligible with respect to the temporal discretization error, and make step time varies. Tables 1 and 2 report the computational results at two selected points t = 1 and t = 2, respectively. We can clearly observe that the order of convergence of temporal discretization is almost O(τ ), which seems reasonable since the employed method in time discretization is of first order.  Table 1. Convergence rates at t = 1 for Example 1 with different step time.
Step time dt The exact solution to this problem is given as u * (x, t) = x 2 /(t + 1) 2 . It can be easily verified that 0 u(t, x) C(T ) for all (x, t) ∈ (−1, 1) × (0, T ) with T < ∞, therefore, the nonlinear source function F is global Lipschitz continuity. As in Example 1, the effectiveness and accuracy should be tested, so we discretize this problem using the approach developed in previous sections using a suitable choice of the parameters discretization M and τ . Figure 2 displays the profiles of the absolute error with different step time τ . As we can see, the approximate and exact solutions match well, which demonstrates that our numerical schema approximates the exact solution greatly.
Next, we check the convergence rates of the temporal discretization. As it done in the previous example, we fix the polynomial degree M = 4, and let step time τ vary as 10 −k . Tables 4 and 5 show that error decreases as the time step τ becomes smaller and as expected convergence order of temporal discretization is O(τ ). That is not surprising since the backward Euler schema used in time discretization is only of first order.  Table 3. Convergence rates at t = 1 for Example 2 with different step time.
Step time dt L 2 -error order L ∞ -error order 10  Table 4. Convergence rates at t = 1.8 for Example 2 with different step time.
Step time dt L 2 -error order L ∞ -error order 10 −1  Table 5. Convergence rates at t = 3 for Example 2 with different step time.
Step time dt L 2 -error order L ∞ -error