Numerical solutions for optimal control problem governed by elliptic system on Lipschitz domains

ABSTRACT

In this paper, an optimal boundary control problem for a distributed elliptic system on Lipschitz domains with boundary homogeneous Dirichlet conditions and independently with Neumann conditions is analysed. The necessary and sufficient optimality conditions for such problems with the quadratic cost functionals are obtained. A Jacobi spectral Galerkin method is introduced to develop a direct solution technique for the numerical solution of elliptic problems subject to Dirichlet and Neumann conditions in one and two dimensions. The numerical examples are included to demonstrate the validity and applicability of the techniques and comparison is made with the existing results. The method is easy to implement and yields very accurate results.

Keywords:

• Optimal control problem
• elliptic system
• Dirichlet and Neumann boundary conditions
• Sobolev spaces on Lipschitz domains
• Jacobi polynomials
• spectral Galerkin method

MSC:

• 46C05
• 49J20
• 93C20

1. Introduction

The boundary value problems on Lipschitz domains have attracted attention in both pure and applied mathematics (see for example [1–3]). On the other hand, there are a lot of engineering problems such as point sensor placements at vertex points of the boundary in the mechanical engineering, corrosive engineering and spacecraft, point controllers at vertex points of the boundary in the stabilization of structural dynamics, etc., which require very careful local analysis around vertex points of the boundary (see [4] and the references therein). The advantage of the Sobolev spaces on Lipschitz domains is connected with the fact that the study of the elliptic partial differential equations becomes very simple and elegant within these spaces. Thus, Sobolev spaces on Lipschitz domains play a very important role in those studies. Most properties of Sobolev spaces on Lipschitz domains are rigorously proved (see e.g. [1,2]).

Optimal control problem for a distributed elliptic system has been a research area of its own (see for example [5–8]). The optimization problem associated with the optimal control of second-order distributed parameter systems defined on smooth domains of ${\mathbb{R}}^n$ has been studied by Lions [5] and recently in [9]. For other cases, we refer for example to [10–12].

Recently, spectral methods have received considerable attention and become increasingly popular for solving different kinds of differential equations. This is definitively due to their superior accuracy for approximating the solutions of problems whose analytical solutions are smooth functions. In practice, this means that highly accurate solutions can be achieved with a small number of discretization nodes. Compared with algebraic convergence rates of the more traditional, finite element, volume, and difference methods spectral methods converge exponentially very fast [13–17].

The spectral methods are divided into Galerkin method, tau method, and collocation method [18–20]. For the convenience of different conditions, three classes of functions are always chosen as the basis for the spectral method: a power basis, a polynomial basis and a trigonometric basis. The important characteristic of the Galerkin method is that the expansion solution must individually satisfy the boundary conditions, whereas in collocation and tau methods, they are not necessary. Physical problems almost always be constrained by known boundary conditions which can be fully incorporated in a Galerkin method [21–23]. Finding a fast and precise numerical solutions for elliptic equations are often a demand in the process of solving problems of fluid dynamics [13,24]. The Galerkin method provides the most suitable numerical approximation for such problems. In his context, Doha and Bhrawy [25–27] introduced a new bases of some orthogonal polynomials for solving the two-dimensional Helmholtz and biharmonic problems by the diagonalization process.

In the present article, we consider a distributed elliptic system of second-order partial differential equation on an open set Ω with a Lipschitz boundary Γ in ${\mathbb{R}}^n$. The existence and uniqueness theorems for the Dirichlet and Neumann problems are proved. The optimal control is characterized by the adjoint problem, and then particular properties of optimal control are also obtained. Furthermore, we aim to develop some efficient spectral algorithms based on the Jacobi Galerkin method (JGM) to approximate the linear second-order elliptic problems such that they can be implemented efficiently. Finally, the accuracy of the proposed algorithms is demonstrated by test problems.

The paper is organized as follows. In Section 2, we recall some basic notation and definitions of the Sobolev space on Lipschitz domains. Section 3 is devoted to formulate the Dirichlet and Neumann problems for elliptic systems on Lipschitz domains. Consequently, the distributed control problems for the Dirichlet problem are studied. We also give the necessary and sufficient conditions for the control to be an optimal for the two problems. Finally, we study the boundary control problem for the Neumann problem. In Section 4, we introduce some algorithms for solving second-order elliptic problems in one and two dimensions by using Jacobi Galerkin spectral methods. Also, the proposed methods are applied to several examples.

2. Sobolev spaces on Lipschitz domains

We first give brief discussions on the Sobolev spaces on the whole domains of ${\mathbb{R}}^n$. As usual ${\mathbb{R}}^n$ denotes the n-dimensional real Euclidian space of points $x=(x_1,\dots ,x_n)$. By a multi-index α, we understand a vector $\alpha =({\alpha }_1,\dots ,{\alpha }_n)$ of non-negative integer components ${\alpha }_i$. We set $|\alpha |=\sum _{i=1}^n{\alpha }_i$.

Let Ω be a non-empty subset of ${\mathbb{R}}^n$. If $f\left(x\right)$ is a sufficiently smooth function on Ω, the partial derivative of $f\left(x\right)$ is defined by $D^{\alpha }f\left(x\right)=\frac{{\mathrm{\partial }}^{|\alpha |}f\left(x\right)}{\mathrm{\partial }{x}_1^{{\alpha }_1}\dots .\mathrm{\partial }{x}_n^{{\alpha }_n}},x^{\alpha }=x_1^{{\alpha }_1}\dots x_n^{{\alpha }_n}.$ In the following, we shall recall some basic function spaces:

For any $r\ge 0$, we set $\begin{array}{rl}{C}^r\left(\mathrm{\Omega }\right)& =\{u:D^{\alpha }u\mathrm{exists}\mathrm{and}\mathrm{continuous}\mathrm{on}\mathrm{\Omega }\\ & \mathrm{for}|\alpha |\le r\},\end{array}$ and put $C^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)=\bigcap _{r\ge 0}{C}^r\left(\mathrm{\Omega }\right)$ the space of infinitely differentiable functions on Ω. $C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ denotes the space of all infinitely differentiable functions with a compact support in Ω, where $supp\left(f\right)=\overline{\{x\in \mathrm{\Omega }:f(x)\ne 0\}}$.

Obviously, $C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ is a real vector space and can be turned into a topological vector space by a proper topology. The space $C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ equipped with the following topology is denoted by $\mathcal{D}\left(\mathrm{\Omega }\right)$: a sequence of functions $\left\{f_k\right\}\subset C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ is said to be convergent to a function $f\in C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ in the space $\mathcal{D}\left(\mathrm{\Omega }\right)$ if

1. there is a compact set $K\subset \mathrm{\Omega }$ such that for all k, $supp(f_k-f)\subset K$ and

2. for each multi-index α, we have $\underset{k\to \mathrm{\infty }}{lim}\parallel D^{\alpha }(f_k-f){\parallel }_{\mathrm{\infty }}$=0.

The space, denoted by ${\mathcal{D}}^{\prime }\left(\mathrm{\Omega }\right)$, of all continuous linear functionals on $\mathcal{D}\left(\mathrm{\Omega }\right)$ is called the (Schwarz) distribution space. $L_{loc}^1\left(\mathrm{\Omega }\right)$ denotes the space of all locally integrable functions (defined almost everywhere on Ω). $L_{loc}^1\left(\mathrm{\Omega }\right)$ can be embedded into ${\mathcal{D}}^{\prime }\left(\mathrm{\Omega }\right)$ using the integration as the duality action. Corresponding to every function $f\in L_{loc}^1\left(\mathrm{\Omega }\right)$ there is a distribution $T_f\in {\mathcal{D}}^{\prime }\left(\mathrm{\Omega }\right)$ defined by $T_f\left(\phi \right)={\int }_{\mathrm{\Omega }}f\phi \mathrm{d}x,\phi \in \mathcal{D}\left(\mathrm{\Omega }\right).$ Evidently $T_f$, thus defined, is a linear functional on $\mathcal{D}\left(\mathrm{\Omega }\right)$.

If $f\in C^{|\alpha |}\left(\mathrm{\Omega }\right)$, then integrating by parts $|\alpha |$-times implies that ${\int }_{\mathrm{\Omega }}{D}^{\alpha }f\left(x\right)\phi \left(x\right)=(-1{)}^{|\alpha |}{\int }_{\mathrm{\Omega }}f(x\left)D^{\alpha }\phi \right(x)\mathrm{for}\mathrm{all}\phi \in \mathcal{D}(\mathrm{\Omega }),$

there are no boundary terms, this is because φ has a compact support in Ω and thus φ and its derivatives vanish near $\mathrm{\partial \Omega }$. This identity motivates the following definition of the derivative $D^{\alpha }T$ of a distribution $T\in {\mathcal{D}}^{\prime }\left(\mathrm{\Omega }\right)$: $\left(D^{\alpha }T\right)\phi =(-1{)}^{|\alpha |}T(D^{\alpha }\phi )\mathrm{for}\mathrm{all}\phi \in \mathcal{D}(\mathrm{\Omega }).$ We are now able to define the concept of a function being the weak derivative of another function. Let α be a multi-index and let $u,v\in L_{loc}^1\left(\mathrm{\Omega }\right)$ such that ${\int }_{\mathrm{\Omega }}v\left(x\right)\phi \left(x\right)=(-1{)}^{|\alpha |}{\int }_{\mathrm{\Omega }}u(x\left)D^{\alpha }\phi \right(x)\mathrm{for}\mathrm{every}\phi \in \mathcal{D}(\mathrm{\Omega }).$

Then the function v is called the weak or distributional α th derivative of u and is denoted by $D^{\alpha }u\left(x\right)$. i.e. $v\left(x\right)=D^{\alpha }u\left(x\right)$.

With the notion of distributions, it is easy now to give a definition for the Sobolev spaces in a form needed in our purpose.

Definition 2.1

Let $m\ge 0$ be an integer. Denote by $H^m\left(\mathrm{\Omega }\right)$, the Sobolev space, the space of all distributions u defined on Ω such that $D^{\alpha }u\in L^2\left(\mathrm{\Omega }\right)\mathrm{for}\mathrm{all}|\alpha |\le m.$ with a finite norm $\parallel u{\parallel }_{m,\mathrm{\Omega }}$ given by $\parallel u{\parallel }_{m,\mathrm{\Omega }}=\sum _{|\alpha |\le m}\parallel D^{\alpha }u{\parallel }_0^2,$ where $||\cdot ||$ is the usual $L^2$-norm.

Lemma 2.1

The Sobolev space $H^m\left(\mathrm{\Omega }\right)$ a Hilbert space together with the inner product $(u,v{)}_m=\sum _{|\alpha |\le m}(D^{\alpha }u,D^{\alpha }v),$ where $(\cdot ,\cdot )$ is the scalar product in $L^2\left(\mathrm{\Omega }\right)$.

Negative order Sobolev spaces: The closure of $C_0^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ with respect to the topology induced by the norm $\parallel \cdot {\parallel }_{m,\mathrm{\Omega }}$ is denoted by $H_0^m\left(\mathrm{\Omega }\right)$. Roughly speaking for $u\in H_0^1\left(\mathrm{\Omega }\right)$ such that $\phi {|}_{\mathrm{\partial \Omega }}=0$ in an appropriate sense. Except m=0 or $\mathrm{\Omega }={\mathbb{R}}^n$, $H_0^m\left(\mathrm{\Omega }\right)$ is a proper subspace of $H^m\left(\mathrm{\Omega }\right)$. For $m\in \mathbb{N}$, $H^{-m}\left(\mathrm{\Omega }\right)$ is defined as the dual space of $H_0^m\left(\mathrm{\Omega }\right)$. In particular, $H^{-1}\left(\mathrm{\Omega }\right):=\left(H_0^1\right(\mathrm{\Omega }){)}^{\ast }$, and for $u\in H^{-1}\left(\mathrm{\Omega }\right)$ $\parallel u{\parallel }_{-1,\mathrm{\Omega }}=\underset{0\ne v\in H_0^1\left(\mathrm{\Omega }\right)}{sup}\frac{|⟨u,v{⟩}_{H^{-1}\left(\mathrm{\Omega }\right)\times H_0^1\left(\mathrm{\Omega }\right)}|}{\parallel v{\parallel }_{1,\mathrm{\Omega }}},$ where the duality pairing on $H^{-1}\left(\mathrm{\Omega }\right)\times H_0^1\left(\mathrm{\Omega }\right)$ is defined by the formula $⟨\phi ,\psi {⟩}_{H^{-1}\left(\mathrm{\Omega }\right)\times H_0^1\left(\mathrm{\Omega }\right)}={\int }_{\mathrm{\Omega }}\psi \left(x\right)D\phi \left(x\right)\mathrm{d}x$ for $\psi \in H^{-1}\left(\mathrm{\Omega }\right)$ and $\phi \in H_0^1\left(\mathrm{\Omega }\right)$.

Proposition 2.2

Embedding of dual Sobolev spaces

For any bounded domain Ω in ${\mathbb{R}}^n,$ one has the following chain of embeddings: $H_0^1\left(\mathrm{\Omega }\right)\hookrightarrow L^2\left(\mathrm{\Omega }\right)=\left(H^1\right(\mathrm{\Omega }){)}^{\ast }\hookrightarrow H^{-1}(\mathrm{\Omega }).$

Proof.

The first embedding is trivial (by definition). The equality follows from the Riesz representation theorem and $L^2\left(\mathrm{\Omega }\right)$ being a Hilbert space. For the last embedding observe that $\{v:\parallel v{\parallel }_{1,\mathrm{\Omega }}\le 1\}\subset \{v:\parallel v\parallel \le 1\}.$ Thus, $f\in L^2\left(\mathrm{\Omega }\right)$ implies $\underset{\{v:\parallel v\parallel \le 1\}}{sup}(f,v)<\mathrm{\infty },$ which implies $\underset{\{v:\parallel v{\parallel }_{1,\mathrm{\Omega }}\le 1\}}{sup}(f,v)<\mathrm{\infty }.$ The proposition is proved.

Definition 2.2

A bounded domain $\mathrm{\Omega }\subset {\mathbb{R}}^n$ is called Lipschitz if near every boundary point $p\in \mathrm{\partial \Omega }$ there exist a rectangular coordinate system $(x,t)$ ($x\in {\mathbb{R}}^n,t\in \mathbb{R}$) and an open neighbourhood U of p such that $\mathrm{\Omega }\cap U=\left\{\right(x,t)\in U|t>\chi (x\left)\right\}$ for some Lipschitz function χ. i.e. $\chi :{\mathbb{R}}^{n-1}\to \mathbb{R}$ and satisfying the Lipschitz condition (1) $|\chi \left(x\right)-\chi \left(x^{\prime }\right)|\le M|x-x^{\prime }|\mathrm{for}\mathrm{all}x,x^{\prime }\in {\mathbb{R}}^{n-1}.$(1)

The smallest M in which (2.1) holds is called the bound of the Lipschitz constant. By choosing a finite open covering $\left\{U_i\right\}$ of $\mathrm{\partial \Omega }$, the Lipschitz constant for a Lipschitz domain is the smallest M such that the Lipschitz constant is bounded by M in every $U_i$. We recall also that a Lipschitz function is almost everywhere differentiable (see e.g. [11]).

Example 2.3

All the smooth domains are Lipschitz. In particular, a domain with smooth $C^1$-boundary is Lipschitz. A very significant non-smooth domain which is Lipschitz is that every polygonal domain in ${\mathbb{R}}^2$ or polyhedron in ${\mathbb{R}}^3$. A more interesting example is that every convex domain in ${\mathbb{R}}^n$ is Lipschitz. A simple example of non-Lipschitz domain is two polygons touching at one vertex only.

3. Elliptic boundary value problems on Lipschitz domains

3.1. Dirichlet problem for the elliptic system

Let Ω be an open set in ${\mathbb{R}}^n$ with a Lipschitz boundary Γ. Consider the following distributed parameter elliptic system with the homogeneous Dirichlet boundary condition: (2) $y-\mathrm{\Delta }y\left(x\right)=f,x\in \mathrm{\Omega },$(2) (3) $y\in H_0^1\left(\mathrm{\Omega }\right),$(3) where $\mathrm{\Delta }=\sum _{i=1}^n\left(2/x_i^2\right)$ is a second-order self-adjoint elliptic partial differential operator maps $H_0^1\left(\mathrm{\Omega }\right)$ onto $H^{-1}\left(\mathrm{\Omega }\right)$. We associate with Δ the bilinear: $\pi (y,\phi )=(y-\mathrm{\Delta }y,\phi ),y,\phi \in H_0^1\left(\mathrm{\Omega }\right),$ and then (4) $\begin{array}{rl}\pi (y,\phi )& =(y-\mathrm{\Delta }y,\phi )\\ & =\sum _{i,j=1}^n{\int }_{\mathrm{\Omega }}\frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_j}\frac{\mathrm{\partial }\phi \left(x\right)}{\mathrm{\partial }{x}_i}\mathrm{d}x+{\int }_{\mathrm{\Omega }}y\left(x\right)\phi \left(x\right)\mathrm{d}x.\end{array}$(4)

Lemma 3.1

The bilinear form $\pi (y,\cdot )$ is coercive on $H_0^1\left(\mathrm{\Omega }\right)\times H_0^1\left(\mathrm{\Omega }\right)$. i.e. there exists a positive constant C such that (5) $\pi (y,y)\ge C||y|{|}_{1,\mathrm{\Omega }}^{2}forally\in H_0^1\left(\mathrm{\Omega }\right).$(5)

Proof.

It is well known that the ellipticity of Δ is sufficient for the coerciveness of $\pi (y,\phi )$ on $H_0^1\left(\mathrm{\Omega }\right)$. By virtue of (4(4) $\begin{array}{rl}\pi (y,\phi )& =(y-\mathrm{\Delta }y,\phi )\\ & =\sum _{i,j=1}^n{\int }_{\mathrm{\Omega }}\frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_j}\frac{\mathrm{\partial }\phi \left(x\right)}{\mathrm{\partial }{x}_i}\mathrm{d}x+{\int }_{\mathrm{\Omega }}y\left(x\right)\phi \left(x\right)\mathrm{d}x.\end{array}$(4) ) and (5(5) $\pi (y,y)\ge C||y|{|}_{1,\mathrm{\Omega }}^{2}forally\in H_0^1\left(\mathrm{\Omega }\right).$(5) ), we have $\begin{array}{rl}\pi (y,y)& ={\int }_{\mathrm{\Omega }}\sum _{i,j=1}^n\frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_j}\frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_i}\mathrm{d}x+{\int }_{\mathrm{\Omega }}y\left(x\right)y\left(x\right)\mathrm{d}x\\ & \ge C{\int }_{\mathrm{\Omega }}\sum _{i=1}^n{\left(\frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_i}\right)}^2\mathrm{d}x+C{\int }_{\mathrm{\Omega }}\left(y\right(x){)}^2\mathrm{d}x\\ & \ge C||y|{|}_{1,\mathrm{\Omega }}^2,\end{array}$ where C is a positive constant.

One can easily show the symmetric property of $\pi (y,\cdot )$, $\pi (y,\phi )=\pi (\phi ,y)$, for all $y,\phi \in H_0^1\left(\mathrm{\Omega }\right)$.

Under the above assumptions, we have the following existence and uniqueness theorem.

Theorem 3.2

Let the associated bilinear form $\pi (y,\cdot )$ satisfies the coerciveness condition (3.4) on $H_0^1\left(\mathrm{\Omega }\right)$. Then, there exists a unique solution to the Dirichlet problem (2(2) $y-\mathrm{\Delta }y\left(x\right)=f,x\in \mathrm{\Omega },$(2) ), (3(3) $y\in H_0^1\left(\mathrm{\Omega }\right),$(3) ).

Proof.

By the assumptions, the bilinear form $\pi (y,\cdot )$ and the linear functional $F\left(\phi \right)=(f,\phi )$ satisfy the hypotheses of the Lax–Milgram theorem for $H_0^1\left(\mathrm{\Omega }\right)$. There exists thus a unique element $u\in H_0^1\left(\mathrm{\Omega }\right)$ such that $\pi (y,\phi )=(y-\mathrm{\Delta }y,\phi )=F\left(\phi \right)\mathrm{for}\mathrm{all}\phi \in H_0^1\left(\mathrm{\Omega }\right).$ This proves the theorem, since y is the weak solution to the Dirichlet problem (2(2) $y-\mathrm{\Delta }y\left(x\right)=f,x\in \mathrm{\Omega },$(2) ), (3(3) $y\in H_0^1\left(\mathrm{\Omega }\right),$(3) ).

3.2. The control problem and the optimization theorem

For a control $u\in L^2\left(\mathrm{\Omega }\right)$, the state $y\left(u\right)$ of the system (2(2) $y-\mathrm{\Delta }y\left(x\right)=f,x\in \mathrm{\Omega },$(2) ), (3(3) $y\in H_0^1\left(\mathrm{\Omega }\right),$(3) ) is given by (6) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f+u\mathrm{in}\mathrm{\Omega }$(6) (7) $y\left(u\right)=0\mathrm{on}\mathrm{\Gamma },$(7) The observation equation is given by $\begin{array}{rl}z\left(u\right)& =M\frac{\mathrm{\partial }}{\mathrm{\partial }{\nu }_{\mathcal{A}}}y\left(u\right),\\ \frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}& =\sum _{i,j=1}^n{a}_{ij}\left(x\right)\cdot \cos (n,x_i)\cdot \frac{\mathrm{\partial }y\left(x\right)}{\mathrm{\partial }{x}_j}\mathrm{on}\mathrm{\Gamma }.\end{array}$ Here, M is a given function in $L^{\mathrm{\infty }}\left(\mathrm{\Gamma }\right)$, n is the exterior normal to Ω at Γ and $\cos (n,x_i)$ is an ith direction cosine of n.

The cost functional $J\left(v\right)$ is given by (8) $\begin{array}{rl}J\left(v\right)& ={\int }_{\mathrm{\Gamma }}\left(M\frac{\mathrm{\partial }}{\mathrm{\partial }{\nu }_{\mathcal{A}}}y\right(v)-z_d{)}^2\mathrm{d}\mathrm{\Gamma }+(Sv,v)\\ & =||M\frac{\mathrm{\partial }}{\mathrm{\partial }{\nu }_{\mathcal{A}}}y\left(v\right)-z_d|{|}_{L^2\left(\mathrm{\Gamma }\right)}^2+(Sv,v),\end{array}$(8) where $z_d$ is a given element in $L^2\left(\mathrm{\Gamma }\right)$ and $S:L^2\left(\mathrm{\Omega }\right)\to L^2\left(\mathrm{\Omega }\right)$ is a linear bounded self-adjoint operator (the solution operator for the problem (6(6) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f+u\mathrm{in}\mathrm{\Omega }$(6) ), (7(7) $y\left(u\right)=0\mathrm{on}\mathrm{\Gamma },$(7) )) satisfying $(Su,u{)}_{L^2\left(\mathrm{\Gamma }\right)}\ge c||u|{|}^2,c>0.$ If ${\mathcal{U}}_{ad}$ (set of admissible controls) is a closed convex subset of $L^2\left(\mathrm{\Omega }\right)$, then making use of Lions's techniques we shall derive the necessary and sufficient conditions of minimizing the cost functional (8(8) $\begin{array}{rl}J\left(v\right)& ={\int }_{\mathrm{\Gamma }}\left(M\frac{\mathrm{\partial }}{\mathrm{\partial }{\nu }_{\mathcal{A}}}y\right(v)-z_d{)}^2\mathrm{d}\mathrm{\Gamma }+(Sv,v)\\ & =||M\frac{\mathrm{\partial }}{\mathrm{\partial }{\nu }_{\mathcal{A}}}y\left(v\right)-z_d|{|}_{L^2\left(\mathrm{\Gamma }\right)}^2+(Sv,v),\end{array}$(8) ). The solving of the formulated optimal control problem is equivalent to seeking a $u\in {\mathcal{U}}_{ad}$ such that (9) $J\left(u\right)=\underset{v\in {\mathcal{U}}_{ad}}{inf}J\left(v\right).$(9)

Theorem 3.3

Under the above assumptions, the problem (9(9) $J\left(u\right)=\underset{v\in {\mathcal{U}}_{ad}}{inf}J\left(v\right).$(9) ) admits a unique solution given by (6(6) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f+u\mathrm{in}\mathrm{\Omega }$(6) ), (7(7) $y\left(u\right)=0\mathrm{on}\mathrm{\Gamma },$(7) ) if, and only if, ${\int }_{\mathrm{\Omega }}\left(p\right(u)+Su)(v-u)\mathrm{d}x\ge 0foreachv\in {\mathcal{U}}_{ad},$ where $p\left(u\right)$ is the adjoint state defined as the solution of the adjoint equation: (10) $p\left(u\right)-\mathrm{\Delta }p\left(u\right)=0in\mathrm{\Omega },$(10) (11) $p\left(u\right)=-M(M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-z_d)on\mathrm{\Gamma }.$(11)

Proof.

It follows from Theorem 3.1 of Lions [5] that the control $u\in {\mathcal{U}}_{ad}$ is optimal if and only if $J^{\mathrm{\prime }}\left(u\right)(v-u)\ge 0\mathrm{for}\mathrm{all}v\in {\mathcal{U}}_{ad}.$ This condition, when explicitly calculated, gives $\begin{array}{rl}& (M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu }-z_d,M\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }\nu }-M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu }{)}_{L^2\left(\mathrm{\Gamma }\right)}\\ & +(Su,v-u)\ge 0\end{array}$ i.e. (12) $\begin{array}{rl}& {\int }_{\mathrm{\Gamma }}(M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu }-z_d,M\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }\nu }-M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu })\mathrm{d}\mathrm{\Gamma }\\ & +(Su,v-u)\ge 0.\end{array}$(12) Multiplying (10(10) $p\left(u\right)-\mathrm{\Delta }p\left(u\right)=0in\mathrm{\Omega },$(10) ) by $\left(y\right(v)-y(u\left)\right)$ and applying Green's formula, we obtain (13) $\begin{array}{rl}0& ={\int }_{\mathrm{\Omega }}\left(p\right(u)-\mathrm{\Delta }p(u\left)\right)\left(y\right(v)-y(u\left)\right)\mathrm{d}x\\ & ={\int }_{\mathrm{\Gamma }}p\left(u\right)(\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }\nu }-\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu })\mathrm{d}\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}\frac{\mathrm{\partial }p\left(u\right)}{\mathrm{\partial }\nu }\left(y\right(v)\\ & -y\left(u\right))\mathrm{d}\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}p(u\left)\right(\left(y\right(v)-y(u\left)\right)\\ & -\mathrm{\Delta }\left(y\right(v)-y(u\left)\right))\mathrm{d}x.\end{array}$(13) In view of (6(6) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f+u\mathrm{in}\mathrm{\Omega }$(6) ) and (7(7) $y\left(u\right)=0\mathrm{on}\mathrm{\Gamma },$(7) ), we have $y\left(v\right)-y\left(u\right)-\mathrm{\Delta }\left(y\right(v)-y(u\left)\right))=v-u$ and $y\left(u\right)=0$ for each $u\in \mathrm{\Gamma }$. Then, the second boundary integral on the right-hand side of (13(13) $\begin{array}{rl}0& ={\int }_{\mathrm{\Omega }}\left(p\right(u)-\mathrm{\Delta }p(u\left)\right)\left(y\right(v)-y(u\left)\right)\mathrm{d}x\\ & ={\int }_{\mathrm{\Gamma }}p\left(u\right)(\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }\nu }-\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu })\mathrm{d}\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}\frac{\mathrm{\partial }p\left(u\right)}{\mathrm{\partial }\nu }\left(y\right(v)\\ & -y\left(u\right))\mathrm{d}\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}p(u\left)\right(\left(y\right(v)-y(u\left)\right)\\ & -\mathrm{\Delta }\left(y\right(v)-y(u\left)\right))\mathrm{d}x.\end{array}$(13) ) vanishes and then it becomes: (14) $0={\int }_{\mathrm{\Gamma }}p\left(u\right)\left(\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}\right)\mathrm{d}\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}p\left(u\right)(v-u)\mathrm{d}x.$(14) Substituting (11(11) $p\left(u\right)=-M(M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-z_d)on\mathrm{\Gamma }.$(11) ) into (14(14) $0={\int }_{\mathrm{\Gamma }}p\left(u\right)\left(\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}\right)\mathrm{d}\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}p\left(u\right)(v-u)\mathrm{d}x.$(14) ), we get $\begin{array}{rl}& {\int }_{\mathrm{\Gamma }}(M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-z_d)(M\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}-M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}})\mathrm{d}\mathrm{\Gamma }\\ & ={\int }_{\mathrm{\Omega }}p\left(u\right)(v-u)\mathrm{d}x\end{array}$ and hence (12(12) $\begin{array}{rl}& {\int }_{\mathrm{\Gamma }}(M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu }-z_d,M\frac{\mathrm{\partial }y\left(v\right)}{\mathrm{\partial }\nu }-M\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }\nu })\mathrm{d}\mathrm{\Gamma }\\ & +(Su,v-u)\ge 0.\end{array}$(12) ) is equivalent to ${\int }_{\mathrm{\Omega }}p\left(u\right)(v-u)\mathrm{d}x+(Su,v-u)\ge 0$ i.e. ${\int }_{\mathrm{\Omega }}\left(p\right(u)+Su)(v-u)\mathrm{d}x\ge 0,$ which completes the proof.

3.3. Neumann problem for the elliptic system

Let us consider the distributed parameter elliptic system with Neumann boundary condition: (15) $y-\mathrm{\Delta }y=f\mathrm{in}\mathrm{\Omega },$(15) (16) $\frac{\mathrm{\partial }y}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=h\mathrm{on}\mathrm{\Gamma }.$(16) The classical result of Meyers and Serrin (cf., e.g. Theorem 3.17 in [1]) enables one to have the following chain of embeddings: $H_0^1\left(\mathrm{\Omega }\right)\hookrightarrow H^1\left(\mathrm{\Omega }\right)\hookrightarrow L^2\left(\mathrm{\Omega }\right).$ Thus, the coerciveness condition (5(5) $\pi (y,y)\ge C||y|{|}_{1,\mathrm{\Omega }}^{2}forally\in H_0^1\left(\mathrm{\Omega }\right).$(5) ) holds on $H^1\left(\mathrm{\Omega }\right)$. Repeating the same arguments as in Section 3.3 and using the trace theorem for the Sobolev spaces, we can prove the following existence and uniqueness theorem for the Neumann problem.

Theorem 3.4

Under the above assumptions, there exists a unique solution $y\in H^1\left(\mathrm{\Omega }\right)$ to the Neumann problem (15(15) $y-\mathrm{\Delta }y=f\mathrm{in}\mathrm{\Omega },$(15) ), (16(16) $\frac{\mathrm{\partial }y}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=h\mathrm{on}\mathrm{\Gamma }.$(16) ).

3.4. Optimization theorem for the boundary control problem

We consider the space $\mathcal{U}=L^2\left(\mathrm{\Gamma }\right)$ (the space of controls), for every control $u\in \mathcal{U}$, the state $y\left(u\right)$ of the system is given as the solution of equation $\pi \left(y\right(u),\varphi )=(u,\varphi )+(h+u,\varphi {)}_{L^2\left(\mathrm{\Gamma }\right),}$ which may be interpreted as (17) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f\mathrm{in}\mathrm{\Omega },$(17) (18) $\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=h+u\mathrm{on}\mathrm{\Gamma },$(18) and the observation is given by $z\left(u\right)=y\left(u\right).$ Finally, the cost function is given by $=||y\left(v\right)-z_d|{|}_{L^2\left(\mathrm{\Omega }\right)}^2+(Sv,v{)}_{L^2\left(\mathrm{\Gamma }\right),}$ where $z_d$ is a given element in $L^2\left(\mathrm{\Omega }\right)$ and $S:L^2\left(\mathrm{\Gamma }\right)\to L^2\left(\mathrm{\Gamma }\right)$ is a linear, bounded, self-adjoint operator satisfying $(Su,u{)}_{L^2\left(\mathrm{\Gamma }\right)}\ge c||u|{|}^2,c>0.$ Now, we have the following optimal control problem: (19) $\underset{v\in {\mathcal{U}}_{ad}}{inf}J\left(v\right),$(19) where ${\mathcal{U}}_{ad}$ (set of admissible controls) is closed convex subset of $\mathcal{U}$. Repeating arguments in the proof of Theorem 3.3, we have the following theorem.

Theorem 3.5

Under the above hypotheses, the optimal control u of (19(19) $\underset{v\in {\mathcal{U}}_{ad}}{inf}J\left(v\right),$(19) ) is characterized by (17(17) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f\mathrm{in}\mathrm{\Omega },$(17) ), (18(18) $\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=h+u\mathrm{on}\mathrm{\Gamma },$(18) ) together with (20) $p\left(u\right)-\mathrm{\Delta }p\left(u\right)=y\left(u\right)-z_{d}in\mathrm{\Omega },$(20) (21) $\frac{\mathrm{\partial }p\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=0on\mathrm{\Gamma },$(21) (22) $\left(p\right(u)+Su,v-u{)}_{L^2\left(\mathrm{\Gamma }\right)}\ge 0forallv\in {\mathcal{U}}_{ad}$(22) where $p\left(u\right)$ is the adjoint state.

Example 3.6

In the case of no constraints on the control (${\mathcal{U}}_{ad}=\mathcal{U}$), we obtain (23) $u=-S^{-1}\left(P\right(u\left){|}_{\mathrm{\Gamma }}\right).$(23) Then, the optimal control is determined by simultaneously solving (17(17) $y\left(u\right)-\mathrm{\Delta }y\left(u\right)=f\mathrm{in}\mathrm{\Omega },$(17) ),(18(18) $\frac{\mathrm{\partial }y\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=h+u\mathrm{on}\mathrm{\Gamma },$(18) ), (20(20) $p\left(u\right)-\mathrm{\Delta }p\left(u\right)=y\left(u\right)-z_{d}in\mathrm{\Omega },$(20) ), (21(21) $\frac{\mathrm{\partial }p\left(u\right)}{\mathrm{\partial }{\nu }_{\mathcal{A}}}=0on\mathrm{\Gamma },$(21) ) (where we eliminate u with the aid of (23(23) $u=-S^{-1}\left(P\right(u\left){|}_{\mathrm{\Gamma }}\right).$(23) )) and then utilizing (23(23) $u=-S^{-1}\left(P\right(u\left){|}_{\mathrm{\Gamma }}\right).$(23) ).

Example 3.7

Consider ${\mathcal{U}}_{ad}=\{u|u\ge 0\mathrm{on}\mathrm{\Gamma }\}$ and $S=\nu \times \mathrm{Id}$, the inequality (22(22) $\left(p\right(u)+Su,v-u{)}_{L^2\left(\mathrm{\Gamma }\right)}\ge 0forallv\in {\mathcal{U}}_{ad}$(22) ) implies $u\ge 0,p\left(u\right)+\nu u\ge 0,u\left(p\right(u)+\nu u)=0\mathrm{on}\mathrm{\Gamma }.$

4. Numerical results and comparisons

We are interested in using the JGM to solve the second-order elliptic equation (24) ${\beta }_{1}u-\mathrm{\Delta }u=f\mathrm{in}\mathrm{\Omega }=I^d,u{|}_{\mathrm{\partial \Omega }}=0,$(24) where $I=(-1,1)$, d=1,2, ${\beta }_1$ is a nonnegative constant, and f is a given source function.

The Jacobi polynomials associated with the real parameters $(\alpha >-1,\beta >-1)$ (see, [28–30]) are a sequence of polynomials $P_n^{(\alpha ,\beta )}\left(x\right)(n=0,1,2,\dots )$, each respectively of degree n, satisfying the orthogonality relation (25) $\begin{array}{rl}& {\int }_{-1}^1(1-x{)}^{\alpha }(1+x{)}^{\beta }{P}_m^{(\alpha ,\beta )}\left(x\right)P_n^{(\alpha ,\beta )}\left(x\right)\mathrm{d}x\\ & =\left\{\begin{array}{ll}0,& m\ne n,\\ h_n,& m=n,\end{array}\right.\end{array}$(25) where $h_n=\frac{2^{\alpha +\beta +1}\mathrm{\Gamma }(n+\alpha +1)\mathrm{\Gamma }(n+\beta +1)}{(2n+\alpha +\beta +1)n!\mathrm{\Gamma }(n+\alpha +\beta +1)}.$ These polynomials are eigenfunctions of the following singular Sturm–Liouville equation: $\begin{array}{rl}& (1-x^2){\phi }^{\mathrm{\prime \prime }}\left(x\right)+[\beta -\alpha -(\alpha +\beta +2\left)x\right]{\phi }^{\mathrm{\prime }}\left(x\right)\\ & +n(n+\alpha +\beta +1)\phi \left(x\right)=0.\end{array}$ Let us first introduce some basic notation that will be used in the sequel. We set $\begin{array}{rl}{S}_N& =\mathrm{span}\left\{P_0^{(\alpha ,\beta )}\right(x),P_1^{(\alpha ,\beta )}(x),\dots ,P_N^{(\alpha ,\beta )}(x\left)\right\},\\ V_N& =\{\nu \in S_N:\nu (\pm 1)=0\},\end{array}$ then the standard Jacobi–Galerkin approximation to (24(24) ${\beta }_{1}u-\mathrm{\Delta }u=f\mathrm{in}\mathrm{\Omega }=I^d,u{|}_{\mathrm{\partial \Omega }}=0,$(24) ) is to find $u_N\in V_N^d$ such that (26) ${\beta }_1(u_N,\nu {)}_w-(\mathrm{\Delta }{u}_N,\nu {)}_w=(f,\nu {)}_w\mathrm{\forall }\nu \in V_N^d,$(26) where $\omega \left(\mathrm{x}\right)=\prod _{i=1}^d(1-x_i{)}^{\alpha }(1+x_i{)}^{\beta }$ and $(u,\nu {)}_w={\int }_{\mathrm{\Omega }}u\nu \omega d\mathrm{x}$ is the scalar product in the weighted space $L_{\omega }^2\left(\mathrm{\Omega }\right)$. The norm in $L_{\omega }^2\left(\mathrm{\Omega }\right)$ will be denoted by $\parallel .{\parallel }_{\omega }$.

Let us denote $H_{\omega }^s\left(\mathrm{\Omega }\right)$ as the weighted Sobolev spaces with the norm $\parallel \cdot {\parallel }_{s,\mathrm{\Omega },\omega }$. It is well known (see [13]) that for ${\beta }_1\ge 0,s\ge 1$, and $u\in H_{\omega }^s\left(\mathrm{\Omega }\right)$, the following optimal error estimate holds: (27) $\parallel u-u_N{\parallel }_{\omega }+N\parallel u-u_N{\parallel }_{1,\mathrm{\Omega },\omega }\le C\left(s\right)N^{-s}\parallel u{\parallel }_{s,\mathrm{\Omega },\omega }.$(27)

In practice, the proper choice of a basis for $V_N^d$ is the key to achieve the optimal convergence rate (26(26) ${\beta }_1(u_N,\nu {)}_w-(\mathrm{\Delta }{u}_N,\nu {)}_w=(f,\nu {)}_w\mathrm{\forall }\nu \in V_N^d,$(26) ) and to obtain a simple linear system as possible. In this section, we introduce some numerical results based on some recent works [25,26,31]. We consider the following examples.

Example 4.1

Consider the one-dimensional Neumann problem $\frac{3}{2}u\left(x\right)-u^{″}\left(x\right)=f\left(x\right)inI,$ with exact solution $u\left(x\right)={\mathrm{e}}^x.$

Similar problems were also investigated by Auteri et al. [18] using a Legendre–Galerkin method (LGM [18]). Table 1 lists the $L^2$- and $H^1$-errors, using the JGM with various choices of N. The numerical results for this problem show that the JGM converges exponentially. We contrast our results with the corresponding results for LGM [18] which we have presented in the fifth and seventh columns of this table.

Table 1. One-dimensional problem with Neumann conditions.

N	α	β	$L^2-\mathrm{error}$		$H^1-\mathrm{error}$
			JGM	LGM [18]	JGM	LGM [18]
	0.5	0.5	$3.63\times {10}^{-5}$		$9.35\times {10}^{-5}$
	0.5	−0.5	$8.59\times {10}^{-5}$		$1.40\times {10}^{-4}$
	0.0.1	0.01	$6.68\times {10}^{-6}$		$5.39\times {10}^{-5}$
4	0	0	$6.68\times {10}^{-6}$	$2.22\times {10}^{-4}$	$5.39\times {10}^{-5}$	$1.88\times {10}^{-3}$
	$-0.5$	0.5	$8.88\times {10}^{-5}$		$1.43\times {10}^{-4}$
	$-0.5$	−0.5	$3.58\times {10}^{-5}$		$9.59\times {10}^{-5}$
	$-0.99$	−0.99	$6.23\times {10}^{-5}$		$1.64\times {10}^{-4}$
	0.5	0.5	$6.19\times {10}^{-10}$		$1.47\times {10}^{-9}$
	0.5	−0.5	$1.39\times {10}^{-9}$		$2.18\times {10}^{-9}$
	0.01	0.01	$\mathrm{4.24.}\times {10}^{-11}$		$5.86\times {10}^{-10}$
8	0	0	$4.15\times {10}^{-11}$	$7.91\times {10}^{-9}$	$5.85\times {10}^{-10}$	$1.34\times {10}^{-7}$
	$-0.5$	0.5	$1.40\times {10}^{-9}$		$2.19\times {10}^{-9}$
	$-0.5$	−0.5	$3.54\times {10}^{-10}$		$1.01\times {10}^{-9}$
	$-0.99$	−0.99	$4.38\times {10}^{-10}$		$1.32\times {10}^{-9}$
	0.5	0.5	$1.02\times {10}^{-16}$		$1.89\times {10}^{-16}$
	0.5	−0.5	$3.14\times {10}^{-16}$		$3.46\times {10}^{-16}$
	0.01	0.01	$1.75\times {10}^{-16}$		$1.92\times {10}^{-16}$
16	0	0	$2.41\times {10}^{-16}$	$3.10\times {10}^{-14}$	$3.16\times {10}^{-16}$	$1.25\times {10}^{-13}$
	$-0.5$	0.5	$1.64\times {10}^{-16}$		$1.85\times {10}^{-16}$
	$-0.5$	−0.5	$2.26\times {10}^{-16}$		$2.51\times {10}^{-16}$
	$-0.99$	−0.99	$5.40\times {10}^{-16}$		$6.00\times {10}^{-16}$

Example 4.2

${\gamma }_{1}u(x,y)-\mathrm{\Delta }u(x,y)=f(x,y)in\mathrm{\Omega }=I\times I,u{|}_{\mathrm{\partial \Omega }}=e$ with exact solution $u(x,y)=x^2+{\mathrm{e}}^{2x+y}$.

This problem is solved in [32] using an LGM in which the integrals on the right-hand side in the resulting linear system are approximated using $(N+1)$-point Legendre Gauss quadrature. The maximum pointwise errors of the proposed JGM using a tensor product process are presented in Table 2. We should note that for all values of N, the present method is always more accurate than the result of LGM [32], which shows the spectral accuracy of our method.

Table 2. $L^{\mathrm{\infty }}$-error of  $u-u_N$  for  $N=8,12,16$.

N	α	β	JGM
	0.5	0.5	$5.87\times {10}^{-4}$
16	0	0	$7.64\times {10}^{-4}$
	−0.5	−0.5	$1.05\times {10}^{-3}$
	0.5	0.5	$2.91\times {10}^{-8}$
24	0	0	$3.95\times {10}^{-8}$
	−0.5	−0.5	$5.18\times {10}^{-8}$
	0.5	0.5	$1.01\times {10}^{-13}$
32	0	0	$1.42\times {10}^{-13}$
	−0.5	−0.5	$1.92\times {10}^{-13}$

It should be noted that the Dirichlet boundary conditions in this example are nonhomogeneous. In [32], the nonhomogeneous Dirichlet boundary conditions at the corners are dealt with using collocation. Then, the nonhomogeneous Dirichlet boundary conditions on each side, excluding the endpoints, are treated using the $L^2$ projection onto the space of polynomials vanishing at the endpoints. In [31], the nonhomogeneous Dirichlet boundary conditions are handled by determining two functions $u^1$ and $u^2$, defined on $\mathrm{\partial \Omega }$, such that $u=u^1+u^2$ on $\mathrm{\partial \Omega }$.

Example 4.3

Consider the problem $u(x,y)-\mathrm{\Delta }u(x,y)=f(x,y),\mathrm{in}\mathrm{\Omega }=I\times I,$ with exact solution $u=x^2+{\mathrm{e}}^{x+2y},$ and the four sets of nonhomogeneous boundary conditions:

• Case 1: $\begin{array}{rl}u(1,y)& ={\mathrm{e}}^r\left(y\right),u(-1,y)={\mathrm{e}}^{\ell }\left(y\right),\\ u_y(x,1)& ={\mathrm{e}}^t\left(x\right),u_y(x,-1)={\mathrm{e}}^b\left(x\right),\end{array}$

• Case 2: $\begin{array}{rl}{u}_x(1,y)& ={\mathrm{e}}^r\left(y\right),u_x(-1,y)={\mathrm{e}}^{\ell }\left(y\right),\\ u(x,1)& ={\mathrm{e}}^t\left(x\right),u(x,-1)={\mathrm{e}}^b\left(x\right),\end{array}$

• Case 3: $\begin{array}{rl}{u}_x(1,y)& ={\mathrm{e}}^r\left(y\right),u_x(-1,y)={\mathrm{e}}^{\ell }\left(y\right),\\ u_y(x,1)& ={\mathrm{e}}^t\left(x\right),u_y(x,-1)={\mathrm{e}}^b\left(x\right),\end{array}$

• Case 4: $\begin{array}{rl}u(1,y)+u_x(1,y)& ={\mathrm{e}}^r\left(y\right),\\ 2u(1,y)-u_x(-1,y)& ={\mathrm{e}}^{\ell }\left(y\right),\\ 2u_y(x,1)-u_y(x,1)& ={\mathrm{e}}^t\left(x\right),\\ -u_y(x,-1)+u_y(x,-1)& ={\mathrm{e}}^b\left(x\right),\end{array}$

In Table 3, we present the maximum pointwise errors of  $u-u_N$ using the JGM with various choices of $\alpha ,\beta ,N$ for the four cases of boundary conditions.

Table 3. $L^{\mathrm{\infty }}$-error of  $u-u_N$  for  $N=8,12,16$.

N	α	β	Case 1	Case 2	Case 3	Case 4
	0.75	0.75	$3.66\times {10}^{-6}$	$2.50\times {10}^{-7}$	$7.48\times {10}^{-6}$	$1.49\times {10}^{-5}$
	0.5	0.5	$2.28\times {10}^{-6}$	$2.08\times {10}^{-7}$	$4.67\times {10}^{-6}$	$9.08\times {10}^{-6}$
8	0	0	$4.08\times {10}^{-7}$	$1.88\times {10}^{-7}$	$2.47\times {10}^{-7}$	$2.57\times {10}^{-7}$
	−0.5	−0.5	$1.15\times {10}^{-6}$	$2.90\times {10}^{-7}$	$2.36\times {10}^{-6}$	$5.36\times {10}^{-6}$
	−0.75	−0.75	$1.45\times {10}^{-6}$	$3.70\times {10}^{-7}$	$2.95\times {10}^{-6}$	$6.52\times {10}^{-6}$
	0.75	0.75	$2.00\times {10}^{-10}$	$1.03\times {10}^{-11}$	$4.05\times {10}^{-10}$	$8.49\times {10}^{-10}$
	0.5	0.5	$1.14\times {10}^{-10}$	$7.47\times {10}^{-12}$	$2.32\times {10}^{-10}$	$4.80\times {10}^{-10}$
12	0	0	$1.20\times {10}^{-11}$	$6.20\times {10}^{-12}$	$7.21\times {10}^{-12}$	$7.10\times {10}^{-12}$
	−0.5	−0.5	$4.40\times {10}^{-11}$	$9.10\times {10}^{-12}$	$8.96\times {10}^{-11}$	$2.00\times {10}^{-10}$
	−0.75	−0.75	$4.92\times {10}^{-11}$	$1.11\times {10}^{-11}$	$9.99\times {10}^{-11}$	$2.20\times {10}^{-10}$
	0.75	0.75	$1.10\times {10}^{-14}$	$6.38\times {10}^{-15}$	$9.18\times {10}^{-15}$	$1.64\times {10}^{-14}$
	0.5	0.5	$1.16\times {10}^{-14}$	$6.21\times {10}^{-15}$	$9.16\times {10}^{-15}$	$1.74\times {10}^{-14}$
16	0	0	$9.83\times {10}^{-15}$	$7.10\times {10}^{-15}$	$7.59\times {10}^{-15}$	$1.33\times {10}^{-14}$
	−0.5	−0.5	$1.10\times {10}^{-14}$	$5.32\times {10}^{-15}$	$8.98\times {10}^{-15}$	$2.52\times {10}^{-14}$
	−0.75	−0.75	$1.44\times {10}^{-14}$	$6.62\times {10}^{-15}$	$7.52\times {10}^{-15}$	$5.20\times {10}^{-14}$

Auteri et al. [18] proposed a direct spectral Galerkin-Legendre solver for the Neumann problem (Case 3) based on a double diagonalization process in which the Gauss-Legendre quadrature formula with N+1 points is used to evaluate the integrals defining the perturbed Neumann condition on the four sides and the term involving the integral of f. The corresponding integrals in the proposed method are evaluated exactly using the Jacobi polynomials.

Acknowledgments

The authors give this work to the spirit of Prof. Dr. A. H. Bhrawy for his helping during the work of this paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work was supported by Taibah University [10.13039/501100002403].

ORCID

G. M. Bahaa http://orcid.org/0000-0001-6670-3110

S. Khidr http://orcid.org/0000-0002-4370-8653