Polynomial-based mean weighted residuals methods for elliptic problems with nonlocal boundary conditions in the rectangle

Abstract. In this paper, polynomial-based mean weighted residuals methods for the solution of elliptic problems with nonlocal boundary conditions, in rectangular domains, are numerically studied. When these methods are employed, the nonclassical boundary conditions involve the solution of large systems of linear equations or least squares problems, hence some numerical techniques for these solvers are compared to show the importance of using efficient algorithms for this purpose. Different kinds of nodes are used to demonstrate how they can be employed to solve different numerical problems when large derivatives of the solution appear. We will also study how using extra precision and/or oversampling can often reduce the computational effort. These methods can also be combined with others (as, for example, finite difference or spline methods) to solve linear and nonlinear parabolic or hyperbolic partial differential equations. The numerical study of some techniques as those explained above can help to obtain significantly better numerical approximations with a smaller computational effort.


Introduction
Nonlocal problems are widely used for mathematical modelling of various processes of physics, thermoelasticity, ecology, chemistry, semiconductor modelling, medicine and biotechnology [1][2][3][4][5], when it is impossible to determine the boundary or initial values of the unknown function and also with problems related with control or inverse PDEs (see [6,7]).Nonclassical boundary and initial-boundary value problems with integral and/or discrete nonlocal boundary conditions were studied for numerous equations.In last decades, numerical methods for the solution of partial differential equations (PDEs) with nonlocal conditions have been developed and analyzed very actively (for example, see [8][9][10][11][12] and references therein).Most of the times, the attention is paid to problems with nonlocal boundary (spatial) conditions.There exist only few studies related with the numerical solution of PDEs with nonlocal initial conditions.c Vilnius University, 2014 Finite difference schemes for two-dimensional elliptic problems with nonlocal boundary conditions were considered in [13][14][15] and references therein.In [16,17], the solvability of the finite difference schemes for a similar problem was proved and the discretization error estimates were obtained.Some error estimates on the finite element approximation for two-dimensional elliptic problem with nonlocal boundary conditions were given in [18].A new method for the solution of the two-dimensional Poisson equation with Bitsadze-Samarskii nonlocal boundary condition was proposed in [19].Various types of finite difference schemes for the solution of nonlocal boundary value problems for multidimensional elliptic equations with nonlocal boundary condition was studied in [20].In [21,22], stable finite difference schemes for the solution of multidimensional elliptic equations with multipoint nonlocal boundary conditions were proposed and analysed.In [23], a radial basis function collocation technique for the solution of a multidimensional elliptic equation with classical Dirichlet boundary condition and nonlocal integral conditions is proposed.However, only in this latter article, problems are formulated on irregular domains.
In this paper, we shall consider that, for all (x, t) ∈ Ω, subject to classical Dirichlet conditions u(x e , y) = 0, y 0 y y e , u(x, y 0 ) = 0, u(x, y e ) = 0, and non-classical boundary condition again with In what follows we will assume that a i (x, y) > 0 for all (x, y) ∈ Ω (i = 0, 1, 2), thus also hence the operator above is elliptic.We assume that a i and f are known and sufficiently smooth functions.
Similar problems have been studied in [24][25][26] and references therein.Actually, some theoretical results obtained in [24,27] can be used partially to study the existence and uniqueness of the weak solution of the variational problem associated to (1)-( 4) or ( 1)-( 3) and (5).However, in this work, we will focus in the numerical implementation of methods to solve it accurately.
Note that it will be possible to extend them to three dimensional elliptic equations similar to equations (1.1) and (1.2) in [28].
It is not necessary to have homogenous conditions (equations ( 2)-( 5)).If the boundary conditions are not homogeneous, simple transformations can be employed to obtain the conditions that appear above.
The organization of this paper is as follows.In Section 2, we will detail the notation of the variational problems associated and introduce a collocation method to solve equation (1) with the nonlocal boundary conditions.The numerical accuracy depends on the different possible of collocation points.In Section 3, some numerical results are provided to support these results and discuss about when the numerical solvers can require extra precision to solve the linear equations or the least square problem.Finally, in Section 4, some conclusions and future work are given.

The variational problem
Let us denote and the Sobolev spaces where α = (α 1 , α 2 ), α i being a natural number and |α| = α 1 + α 2 , We use the traditional scalar products and norms: Observe that, in this case, when we try to find the variational formulation of the problem and according to Green's formula the last part does not vanish in Γ 2 .Additionally, in this case, we do not necessarily have W ⊇ H 1 0 (Ω) nor u has to be a function in H 1 0 (Ω) either.The existence and uniqueness of the solution of the variational problem has been considered previously (for example, see [24] and references therein for some cases) without deep numerical studies.Hence, here we focus on its numerical study with the collocation pseudo-spectral method.

The weighted residual method
The collocation method can be considered as a generalized Galerkin method and therefore it goes under the stability and convergence analysis of these methods.
Thus, we will work in the Hilbertian basis of the Cartesian product of polynomials on x and y, i.e., the solution will be where the coefficients a kl s are to be determined.This solution can be approximated at the N th order by the dot product in what follows we will use Nonlinear Anal.Model.Control, 2014, Vol.19,No. 3,[448][449][450][451][452][453][454][455][456][457][458][459] These coefficients are determined imposing that the residual of the discretization, R 1 (x, y, a k,l ) 2 + R 2 (x, y, a k,l ) 2 at certain nodes {(x 0 , y 0 ), (x 0 , y 1 ), . . ., (x M , y M )} is minimum.Note that M and N do not need to be equal.When M = N , methods for solving large systems of linear equations are used.In general, singular value decompositions and iterative methods are considered for these problems, however, many of these systems are also ill-conditioned [29] and therefore preconditioning techniques such as preconditioned conjugate gradient (PCG), m-step Jacobi PCG or the preconditioned GMRES methods are considered.
When N < M , we can consider the collocation method as a least square mean weighted residual method, i.e., the least squares determines the coefficients a k,l in such a way that the inner product of R(x; y; a k,l ) with itself is made as small as possible for the functions in (17).In this work, we will use the LeastSquares command of MATHEMATICA.
As for the grid points we will consider the cartesian product of equispaced and also Chebyshev nodes (x 0 = y 0 = 0, x e = x M = y M = 1 for the boundary conditions).For complex geometries, equispaced points usually are the better choice.However, for rectangles as it is the case, some pseudospectral grids are more efficient (see [30]).Chebyshev nodes minimize the maximal possible error and solve the famous Runge phenomenon.Additionally, the cost of some calculus can be reduced through fast fourier transforms (FFTs) in some problems (however, the reader should remember that in this paper there are also nonlocal boundary conditions, not only traditional homogeneous conditions are considered).

Numerical examples
Example 1.Let us first consider problem (1)-( 3) with nonlocal condition (5), when In this case, the exact solution is u(x, y) = (x − 1)(y − 1) sin(πy/2).We used both, equispaced and Chebyshev nodes without, Γ = 1, and with oversampling, Γ = 3/2, 2 and 5/2.In Fig. 1, we show the errors only when Γ = 1 and 2 (the results for the other two values were very similar).The errors are calculated with the infinity norm (maximum norm) at a uniform and dense grid.We used the LeastSquares solver (a direct method) of MATHEMATICA to solve the least squares problem.
First of all, we see that using these techniques we obtain a very fast convergence, but also that mean weighted residuals methods have some limitations, in general, with standard precision, with equispaced points we are not able to obtain errors smaller than  We can usually obtain better results when using Chebyshev points than with equispaced nodes and it is frequently better to use oversampling, both techniques are specially useful when large values of N are required.Still there are problems where differences are not big or even more, where equispaced nodes achieve slightly better results (and something similar happens between using, or not, oversampling).Why, then, so many researchers used Chebyshev and other kinds of grids when they solved PDEs with classical conditions?First, because Chebyshev nodes are very important in approximation theory.They form a particularly good set of nodes for polynomial interpolation: given a function f on the interval [0, 1] and M + 1 points x 0 , x 1 , . . ., x M , in that interval, there is only a unique polynomial P M of degree M which has value f (x i ) at each point x i and the interpolation error at x is for some ξ in [0, 1].The Chebyshev nodes are the ones which minimize (and therefore are more efficient when N, M → ∞).
Even more, in some cases equispaced points (when Γ = 1) do not always guarantee convergence (a famous example in interpolation in one dimension is the Runge phenomenon, 1/(1 + a 2 x 2 ), with a → ∞), whilst Chebyshev nodes have always better convergence for these kinds of equations with very large derivatives of the function to interpolate.The reason is these points cluster quadratically near the end points 0 and 1. Obviously, other kinds of nodes (as, for example, some Legendre or Gaussian quadratures) guarantee convergence for such kind of problems, but they require to cluster quadratically near the end points.It is also possible to obtain convergence in such interpolation problems when we impose M = O(N 2 ), but this technique can become very expensive then.On the other hand, if the region where we are solving the PDE is complex, frequently equispaced points need to be used.
When the PDE has traditional boundary conditions, another procedure is splitting the complex regions in a set of other simpler regions, for example, in triangles or tetrahedra (and using some quadrature nodes there).This could be quite more challenging when PDEs have nonlocal boundary conditions since (4) or ( 5) would be an equation depending on the coefficients a k,l for several triangles or tetrahedra (still this could be solved using extra conditions, for example, given by some oversampling).Obviously, it could be interesting and much more difficult, for example, the study of the decomposition of complex regions when the PDE has more than one nonstandard boundary condition.
Example 2. For the second numerical example, let us consider u(x, y) = e −x x 2 + y + xy 2 , hence the boundary conditions cannot be homogeneous, actually, the equations now are We already stated that the convergence analysis of these schemes undergoes the properties of Galerkin methods and the theoretical study of these algorithms has not been provided for boundary conditions as those in this test problem.However, collocation schemes have been demonstrated efficient with two-point BVPs and also with some PDEs without homogeneous conditions (see [31,32], for example).This can make the reader thinks that these procedures will work in these kinds of examples.Effectively, they do as we can check in Fig. 3 (the Robin boundary condition does not affect either).There we show the errors obtained with Γ = 1 and Γ = 2 for equispaced and Chebyshev nodes, without and with oversampling.Again errors in standard precision are over O(10 −11 )-O(10 −12 ).
We can solve this inconvenient with extra precision.In the left part of Fig. 4, we check the fast convergence of these methods when the conditions are not homogeneous and how Chebyshev nodes provide better results in general and particularly when N, M → ∞.In the right part of the figure we show the errors when N = 12, Γ = 1 with standard precision and Chebyshev nodes this time.
Example 3.Many papers have been recently published about the resolution of highly oscillatory problems which can become a quite challenging topic ( [33][34][35] can provide surveys and references about different techniques, several kinds of quadratures, radial basis function and wavelets, to solve oscillatory ODEs and PDEs).The third numerical example presents a mildly oscillatory PDE with nonlocal condition.
First of all, we see how oscillatory problems require several nodes per wavelength to resolve oscillations and obtain "small" errors.Additionally oscillations, and therefore large derivatives make that errors, in finite precision, are O(10 −5 ) in this problem.This is enough for many engineering problems, but could be worrying for many researchers, specially because sin(ωy), with ω = 3π in the solution, can be considered a not very highly oscillatory function.We can understand that, if ω 10 (for example), it could difficult to obtain small enough errors.However, in the right-side figure, we can check that using values of N big enough an extra precision we can solve this problem.In this kind of problems, it is clear that using Chebyshev nodes or oversampling can be an important tool to obtain more accurate numerical results.
In this paper, we are not studying the comparison among direct and iterative solvers for the least squares problem.However, the reader can understand from this numerical example that this could be a very important issue for many PDEs.In this paper, we just considered direct methods for two-dimensional PDEs.However, the implementation of such kind of pseudospectral schemes can be done in a similar way for three-dimensional www.mii.lt/NAelliptic PDEs; but, for some kind of problems where large values of N and/or extra precision are required to obtain accurate results (as oscillatory problems), these schemes could become very expensive, hence iterative solvers should be studied.

Conclusions
In this paper, polynomial-based mean weighted residuals methods have been proposed to solve elliptic problems with nonlocal boundary conditions.The paper focuses on the numerically study of such schemes.Different kinds of nodes and oversampling have been used to demonstrate how they can be employed to solve different numerical problems.We also studied how using extra precision is necessary for some kind of problems whenever we want to obtain small enough errors.These techniques can easily be extended to three dimensional PDEs with nonlocal boundary conditions.However, in such cases, the dimensions of the least squares are larger and the study of iterative methods can become even more interesting.The convergence rate of most of these iterative schemes is related with the conditioning of the matrices (related with the linear systems we solve), hence it would be necessary to study the smallest singular values of these matrices and preconditioning techniques (see [30,36]).
Obviously, they can also be used to solve nonlinear equations as it was done in [32] with hyperbolic nonlinear Klein-Gordon equations (in this paper, traditional conditions were considered).Normally Newton-type schemes are utilized in these cases.In [32], moderate values of N were considered only, obtaining good results, however, the convergence of these algorithms when N is larger could be slower.Additionally, if the solutions of the PDEs have large derivatives it could be necessary the use of extra working precision (in a similar way as we show in the third numerical test).Both issues (more working precision and slower convergence) could make that the computational effort clearly grows to solve these nonlinear equations.For nonlinear hyperbolic and parabolic equations, some spectral collocation methods have also been combined with Runge-Kutta methods, finite-difference or spline methods [37,38], but with neither highly oscillatory nor other PDEs with large derivatives therein.This should be deeply analyzed in another paper.
Another possible goal for the future can be the study of similar pseudospectral collocation methods in other kind of geometries as a triangle or a tetrahedron, but also the circle or more irregular domains.The techniques in this paper can be easily used in this kind of geometries.However, in some of these cases, it could be necessary the use of oversampling to obtain accurate numerical results.

Fig. 2 .
Fig. 2. The errors for different values of N for equispaced (squares) and Chebyshev nodes (circles) with Γ = 3/2 with extra precision (left).In the right-side part, a three-dimensional plot of the errors is shown with N = 10, Γ = 3/2 with standard precision and equispaced nodes are used.O(10 −11 )-O(10 −12 ).These limitations are not produced by the nonlocal conditions, they also appear with homogeneous conditions.They are produced, firstly, because the operator uses singular value decompositions with truncated parameter.Additionally, we have the errors at any other point in the interval produced by the M + 1 derivatives of the function f , plus O(M 2 ) errors due to the second derivatives in the PDE.These limitations disappear when extra precision (and extra digits in the tolerance of the LeastSquares command) is used as we can check in the left part of Fig.2.In the right part of the figure, we show the errors when N = 10, Γ = 3/2 with standard precision and equispaced nodes.When we use extra precision, we can clearly check the spectral convergence of these methods.We can usually obtain better results when using Chebyshev points than with equispaced nodes and it is frequently better to use oversampling, both techniques are specially useful when large values of N are required.Still there are problems where differences are not big or even more, where equispaced nodes achieve slightly better results (and something similar happens between using, or not, oversampling).Why, then, so many researchers used Chebyshev and other kinds of grids when they solved PDEs with classical conditions?First, because Chebyshev nodes are very important in approximation theory.They form a particularly good set of nodes for polynomial interpolation: given

Fig. 5 .
Fig. 5.The errors for different values of N for equispaced (squares) and Chebyshev nodes (circles) with Γ = 1 with normal precision (left) and extra precision (right).