Explicit Discrete Solution for Some Optimization Problems and Estimations with Respect to the Exact Solution

1. Introduction

We consider a multidimensional bounded domain $\Omega \subset {\mathbb{R}}^n$ whose regular boundary $\Gamma$ consists of three disjoint portions ${\Gamma }_i$ with $meas\left({\Gamma }_i\right)>0,$ for $i=1,2,3$. We define two stationary heat conduction problems $\left(S\right)$ and $\left(S_{\alpha }\right)$ with mixed boundary conditions which are given by (1)-(2) and (1)-(3) respectively:

$$\begin{array}{ccc}& -\Delta u=gin\Omega ,& \hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill u{|}_{{\Gamma }_1}=b,& -\frac{\partial u}{\partial n}{|}_{{\Gamma }_2}=q,& -\frac{\partial u}{\partial n}{|}_{{\Gamma }_3}=0,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill -\frac{\partial u}{\partial n}{|}_{{\Gamma }_1}=\alpha (u-b),& -\frac{\partial u}{\partial n}{|}_{{\Gamma }_2}=q& -\frac{\partial u}{\partial n}{|}_{{\Gamma }_3}=0,\hfill \end{array}$$

(3)

where g is the internal energy of the system in $\Omega$, $b>0$ the environmental temperature on ${\Gamma }_1$, q is the heat flux on ${\Gamma }_2$ and $\alpha >0$ is the convective heat coefficient on ${\Gamma }_1.$ We assume that: $g\in H=L^2(\Omega )$, $q\in Q=L^2\left({\Gamma }_2\right)$ and $b\in B=H^{1/2}\left({\Gamma }_1\right).$ These problems can be considered as Stefan’s stationary problems [8,19,20,21,22]. Notice that mixed boundary conditions play an important role in several applications, e.g. heat conduction and electric potential problems [12].

The variational formulation of the elliptic problems (1)-(2) and (1)-(3) can be found in [8,9,10,19]. It can be seen that they have unique solution [20]. In general, the solutions of mixed elliptic boundary problems is not so regular [11] but there are cases in which they are regular [4,15,18]. Other theoretical optimization problems in the subject have been done in [13,25].

We define the optimization problems $\left(P_i\right)$ and $\left(P_{i\alpha }\right)$$i=1,2,3$, associated to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$, respectively (see [5,9,16,17,24]).

The distributed optimization problems $\left(P_1\right)$ and $\left(P_{1\alpha }\right)$ on the constant internal energy g are formulated as:

$$\begin{array}{cc}\hfill & \mathrm{find}{g}_{op}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_1\left(g_{op}\right)=\underset{g\in \mathbb{R}}{min}{J}_1\left(g\right)\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & \mathrm{find}{g}_{{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{1\alpha }\left(g_{{\alpha }_{op}}\right)=\underset{g\in \mathbb{R}}{min}{J}_{1\alpha }\left(g\right)\hfill \end{array}$$

(5)

where $J_1:\mathbb{R}\to R_0^{+}$ and $J_{1\alpha }:\mathbb{R}\to R_0^{+}$ are given by

$$J_1\left(g\right)={\displaystyle \frac{1}{2}}{∥u_g-z_d∥}_H^2+{\displaystyle \frac{M_1}{2}}{∥g∥}_H^2,J_{1\alpha }\left(g\right)={\displaystyle \frac{1}{2}}{∥u_{\alpha g}-z_d∥}_H^2+{\displaystyle \frac{M_1}{2}}{∥g∥}_H^2$$

(6)

with $M_1\in {\mathbb{R}}^{+}$ and $z_d\in \mathbb{R}$. For each $g\in \mathbb{R}$, the functions $u_g$ and $u_{\alpha g}$ denote the unique solutions to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ respectively, for data $q\in \mathbb{R}$, $b\in \mathbb{R}$.

The boundary optimization problems $\left(P_2\right)$ and $\left(P_{2\alpha }\right)$ on the constant heat flux q on ${\Gamma }_2$ are defined as:

$$\begin{array}{cc}\hfill & \mathrm{find}{q}_{op}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_2\left(q_{op}\right)=\underset{q\in \mathbb{R}}{min}{J}_2\left(q\right)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \mathrm{find}{q}_{{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{2\alpha }\left(q_{{\alpha }_{op}}\right)=\underset{q\in \mathbb{R}}{min}{J}_{2\alpha }\left(q\right)\hfill \end{array}$$

(8)

where $J_2:\mathbb{R}\to R_0^{+}$ and $J_{2\alpha }:\mathbb{R}\to R_0^{+}$ are given by

$$J_2\left(q\right)={\displaystyle \frac{1}{2}}{∥u_q-z_d∥}_H^2+{\displaystyle \frac{M_2}{2}}{∥q∥}_Q^2,J_{2\alpha }\left(q\right)={\displaystyle \frac{1}{2}}{∥u_{\alpha q}-z_d∥}_H^2+{\displaystyle \frac{M_2}{2}}{∥q∥}_Q^2$$

(9)

with $M_2\in {\mathbb{R}}^{+}$ and $z_d\in \mathbb{R}$. For each $q\in \mathbb{R}$, we denote with $u_q$ and $u_{\alpha q}$ the unique solutions to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ respectively, for data $g\in \mathbb{R}$ and $b\in \mathbb{R}$.

The boundary optimization problems $\left(P_3\right)$ and $\left(P_{3\alpha }\right)$ on the constant temperature b in an external neighborhood of ${\Gamma }_1$ are set as

$$\begin{array}{cc}\hfill & \mathrm{find}{b}_{op}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_3\left(b_{op}\right)=\underset{b\in \mathbb{R}}{min}{J}_3\left(b\right)\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \mathrm{find}{b}_{{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{3\alpha }\left(b_{{\alpha }_{op}}\right)=\underset{b\in \mathbb{R}}{min}{J}_{3\alpha }\left(q\right)\hfill \end{array}$$

(11)

where $J_3:\mathbb{R}\to R_0^{+}$ and $J_{3\alpha }:\mathbb{R}\to R_0^{+}$, given by

$$J_3\left(b\right)={\displaystyle \frac{1}{2}}{∥u_b-z_d∥}_H^2+{\displaystyle \frac{M_3}{2}}{∥b∥}_B^2,J_{3\alpha }\left(b\right)={\displaystyle \frac{1}{2}}{∥u_{\alpha b}-z_d∥}_H^2+{\displaystyle \frac{M_3}{2}}{∥b∥}_B^2$$

(12)

with $M_3\in {\mathbb{R}}^{+}$ and $z_d\in \mathbb{R}$. For every $b\in \mathbb{R}$, the functions $u_b$ and $u_{\alpha b}$ are the unique solutions of the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ respectively, for data $g\in \mathbb{R}$ and $q\in \mathbb{R}$, .

In [6] the explicit solutions to the continuous systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ were obtained, as well as the associated optimization problems $\left(P_i\right)$ and $\left(P_{i\alpha }\right)$ for $i=1,2,3$, in the particular case when the domain is a rectangle. These explicit solutions provide a benchmark for testing the accuracy of numerical methods.

The goal of this paper is to find the discrete explicit solutions to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$, through the finite difference method, in a rectangular domain. In addition, we aim to obtain the discrete explicit solutions to the optimization problems $\left(P_i\right)$ and $\left(P_{i\alpha }\right)$ for $i=1,2,3$. Taking advantage of the exact explicit solutions, we will estimate the order of convergence of the approximate solutions as a function of the step discretization.

It is worth to mention that there are several articles available in the literature that obtain explicit discrete solutions of some optimization problems [1,2,3]. For example, in [14], exact formulas are derived for the solution of an optimal boundary control problem governed by the one dimensional heat equation where the control function measures the distance of the final state from the target. In [26] a finite element approximation is applied for some kind of parabolic optimal control problems with Neumann boundary conditions. Some numerical experiments are carried out setting a rectangular domain.

This paper is organized as follows: in Section 2 we obtain the discrete explicit solution to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ by the finite difference method. In Section 3, we obtain explicit discrete solutions to the discrete distributed optimization problems associated to $\left(P_1\right)$ and $\left(P_{1\alpha }\right)$, respectively where the variable is the internal energy g. In Section 4, we define discrete boundary optimization problems where the variable is the heat flux q, associated to $\left(P_2\right)$ and $\left(P_{2\alpha }\right)$, respectively obtaining the discrete explicit solutions. In the same manner, in Section 5, we derive explicit discrete solutions to the discrete boundary optimal control problems associated to $\left(P_3\right)$ and $\left(P_{3\alpha }\right)$, respectively where the optimization variable is b. In all cases, when the step discretization goes to zero, convergence results are obtained estimating also the order of convergence of the approximate solutions. In Section 6, we carry out some numerical simulations in order to illustrate the convergence theoretical results obtained in the previous sections. Finally, in Section 7, we analyze the order of convergence of the discrete systems associated with $\left(S\right)$ and $\left(S_{\alpha }\right)$ by considering a modified approximation of the Neumann boundary condition on ${\Gamma }_2$, which leads to an improved convergence order.

2. Discrete Systems for $\left(S\right)$ and $\left(S_{\alpha }\right)$

In this section we obtain the discrete explicit solutions to the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ in a particular domain.

Let us consider a rectangular domain in the plane $\Omega =(0,x_0)\times (0,y_0)$ with $x_0>0$ and $y_0>0$. Its boundaries ${\Gamma }_i$ for $i=1,2,3$ are defined by:

$${\Gamma }_1=\{(0,y):y\in (0,y_0]\},{\Gamma }_2=\{(x_0,y):y\in (0,y_0]\}$$

and

$${\Gamma }_3=\{(x,0):x\in [0,x_0]\}\cup \{(x,y_0):x\in [0,x_0]\}.$$

According to [6] the continuous solutions, in $\Omega$, for the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ defined by (1)-(2) and (1)-(3) are given by:

$$\begin{array}{c}u(x,y)=-\frac{1}{2}g{x}^2+(g{x}_0-q)x+b,\forall (x,y)\in \Omega \hfill \\ \\ u_{\alpha }(x,y)=-\frac{1}{2}g{x}^2+(g{x}_0-q)x+\frac{1}{\alpha }(g{x}_0-q)+b,\forall (x,y)\in \Omega .\hfill \end{array}$$

(13)

As consequence of the symmetry of the domain $\Omega$, the solution u and $u_{\alpha }$ of the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ are independent of the variable y, and therefore we will work with one-dimensional problems.

Given $n\in \mathbb{N}$, we define:

$$h=\frac{x_0}{n};x_i=(i-1)h,\mathrm{for}i=1,\dots ,n+1,{\mathrm{u}}_i\approx u(x_i,y)\mathrm{for}i=2,\dots ,n+1,$$

(14)

where, n is the number of subintervals to consider in $[0,x_0]$, h is the constant size of these subintervals, ${\mathrm{u}}_1=b$ and ${\mathrm{u}}_i$ is the approximate value of the function u in $(x_i,y)\in \Omega$ for $i=2,\dots ,n+1$.

We apply the classical finite–difference method to the system $\left(S\right)$ described by equations (1)–(2). Since the boundary condition on ${\Gamma }_1$ prescribes $u(0,y)=b$, we immediately obtain ${\mathrm{u}}_1=b$.

For the interior nodes, we use the classical centered second–order finite difference approximation:

$${\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}(x_i,y)\approx {\displaystyle \frac{u(x_{i+1},y)-2u(x_i,y)+u(x_{i-1},y)}{h^2}},i=2,3,\dots ,n,$$

(15)

and from (1) it follows that

$$-g{h}^2={\mathrm{u}}_{i+1}-2{\mathrm{u}}_i+{\mathrm{u}}_{i-1},i=2,3,\dots ,n.$$

(16)

To incorporate the Neumann boundary condition on ${\Gamma }_2$, we use a backward finite difference for the first derivative:

$${\displaystyle \frac{\partial u}{\partial x}}(x_{n+1},y)\approx {\displaystyle \frac{u(x_{n+1},y)-u(x_n,y)}{h}},$$

(17)

which, using the boundary condition ${\displaystyle \frac{\partial u}{\partial x}}(x_{n+1},y)=q$, leads to

$$-qh={\mathrm{u}}_{n+1}-{\mathrm{u}}_n.$$

(18)

As a consequence, we obtain the discrete linear system $\left(S_h\right)$

$$A\overline{\mathrm{v}}=T,$$

(19)

where $\overline{\mathrm{v}}={\left({\mathrm{u}}_i\right)}_{i=2,\dots ,n+1}\in {\mathbb{R}}^n$ is the vector of unknowns, A is the associated tridiagonal coefficient matrix:

$$A=\left(\begin{array}{cccccccc}-2& 1& 0& \dots & \dots & & 0& \\ 1& -2& 1& 0& \dots & & 0\\ 0& 1& -2& 1& 0& \dots & 0\\ ⋮\\ 0& \dots & 0& 1& -2& 1& 0\\ 0& \dots & \dots & 0& 1& -2& 1\\ 0& \dots & \dots & \dots & 0& -1& 1\end{array}\right)$$

(20)

and $T\in {\mathbb{R}}^n$ is the vector of independent terms:

$$T={\left(-g{h}^2-b,-g{h}^2,\dots ,-g{h}^2,-qh\right)}^t.$$

(21)

The square matrix A is invertible and its inverse matrix is given by

$$A^{-1}=\left(\begin{array}{cccccccc}-1& -1& -1& \dots & \dots & -1& 1& \\ -1& -2& -2& \dots & \dots & -2& 2& \\ -1& -2& -3& \dots & \dots & -3& 3& \\ ⋮& ⋮& ⋮& & & ⋮& ⋮\\ -1& -2& -3& \dots & & -(n-1)& n-1\\ -1& -2& -3& \dots & & -(n-1)& n\end{array}\right)$$

Then the system $\left(S_h\right)$ has a unique solution:

$${\mathrm{u}}_1=b{\mathrm{u}}_i=b+h^{2}g\left((i-1)n-\frac{i(i-1)}{2}\right)-(i-1)hq,i=2,\dots ,n+1.$$

As $n=\frac{x_0}{h}$, it follows that:

$${\mathrm{u}}_i=b+(i-1)h(g{x}_0-q)-h^{2}g\frac{i(i-1)}{2},i=1,\dots ,n+1.$$

Then, the continuous solution $u(x,y)$ of system $\left(S\right)$ can be approximated by the piecewise linear interpolant $u_h(x,y)$ obtained from the nodal values computed by the finite difference scheme. More precisely, we define

$$u_h(x,y)=(g{x}_0-q-hgi)x+h^{2}g\left(\frac{i(i-1)}{2}\right)+b,x\in [x_i,x_{i+1}],y\in [0,y_0],$$

(22)

with $i=1,\dots ,n.$

In the following lemma, we give some bounds for the approximate function $u_h:$

Lemma 1. 

a) 

For every $(x_i,y)$ with $i=1,\dots ,n+1$, $y\in [0,y_0]$, we have that:

i) 

if $g>0$ then $u_h(x_i,y)\le u(x_i,y)$.

ii) 

if $g<0,$ then $u_h(x_i,y)\ge u(x_i,y)$.

b) 

The following bounds hold:

$${\parallel u-u_h\parallel }_H\le C_{1}h,\mathit{and}{\parallel \frac{\partial u}{\partial x}-\frac{\partial u_h}{\partial x}\parallel }_H\le \widetilde{C_1}h,$$

where $C_1=x_0\left|g\right|\sqrt{\frac{2}{15}{x}_0{y}_0}$ and $\widetilde{C_1}=\left|g\right|\sqrt{\frac{1}{3}{x}_0{y}_0}.$

Proof. 

a)

From functions u and $u_h$, we have that

$$u(x_i,y)-u_h(x_i,y)=u(x_i,y)-{\mathrm{u}}_i=\frac{g{h}^2(i-1)}{2},i=1,\dots ,n+1.$$

b)

From the definition of norm in space H and the formulas (13) and (22) for the functions u and $u_h$, respectively, it follows that:

$$\begin{array}{cc}\left|\right|u-u_h{\left|\right|}_H^2\hfill & {\displaystyle =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left(u(x,y)-u_h(x,y)\right)}^{2}dx}\hfill \\ \\ & ={\displaystyle \frac{1}{120}}{y}_0{h}^5{g}^2{n}^3\left({\displaystyle \frac{1}{n^2}}+{\displaystyle \frac{5}{n}}+10\right)\hfill \\ \\ & \le {\displaystyle \frac{2}{15}}{y}_0{h}^5{g}^2{n}^3={\displaystyle \frac{2}{15}}{h}^2{g}^2{x}_0^3{y}_0=C_1^2{h}^2.\hfill \end{array}$$

The norm $\left|\right|\frac{\partial u}{\partial x}-\frac{\partial u_h}{\partial x}{\left|\right|}_H$ can be computed analogously.

□

We apply the classical finite–difference method to the system $\left(S_{\alpha }\right)$ defined by equations (1)–(3) where ${\mathrm{u}}_{\alpha ,i}\approx u_{\alpha }(x_i,y)$ for $i=1,\dots ,n+1$.

The Robin boundary condition on ${\Gamma }_1$ is approximated by a forward finite–difference scheme, namely,

$${\displaystyle \frac{{\mathrm{u}}_{\alpha ,2}-{\mathrm{u}}_{\alpha ,1}}{h}}=\alpha \left({\mathrm{u}}_{\alpha ,1}-b\right).$$

(23)

Moreover, at the interior nodes we use the approximation given in (16), while the Neumann boundary condition at $x_{n+1}$ is discretized according to (17).

Then we obtain a linear system $\left(S_{h\alpha }\right)$:

$$A_{\alpha }{\overline{\mathrm{v}}}_{\alpha }=T_{\alpha }$$

where the vector of unknowns ${\overline{\mathrm{v}}}_{\alpha }\in {\mathbb{R}}^{n+1}$ is given by ${\overline{\mathrm{v}}}_{\alpha }={\left({\mathrm{u}}_{\alpha ,i}\right)}_{i=1,\dots n+1}$, the tridiagonal coefficient matrix $A_{\alpha }$ of order $n+1$, is defined as:

$$A_{\alpha }=\left(\begin{array}{cccccccc}-(1+\alpha h)& 1& 0& \dots & \dots & & 0& \\ 1& -2& 1& \\ ⋮& \dots & \dots & & 1& -2& 1\\ 0& 0& 0& \dots & 0& -1& 1\end{array}\right)$$

(24)

and

$$T_{\alpha }={\left(-\alpha bh,-g{h}^2,\dots ,-g{h}^2,-qh\right)}^t\in {\mathbb{R}}^{n+1}.$$

(25)

It can be seen that the square matrix $A_{\alpha }$ is invertible and its inverse matrix is given by

$$A_{\alpha }^{-1}={\displaystyle \frac{1}{\alpha h}}\left(\begin{array}{cccccccc}-1& -1& -1& \dots & \dots & -1& 1& \\ -1& -(1+\alpha h)& -(1+\alpha h)& \dots & \dots & -(1+\alpha h)& 1+\alpha h& \\ -1& -(1+\alpha h)& -(1+2\alpha h)& \dots & \dots & -(1+2\alpha h)& -(1+2\alpha h)& \\ ⋮& ⋮& ⋮& & & ⋮& ⋮\\ -1& -(1+\alpha h)& -(1+2\alpha h)& \dots & & -(1+(n-1\left)\alpha h\right)& 1+(n-1)\alpha h\\ -1& -(1+\alpha h)& -(1+2\alpha h)& \dots & & -(1+(n-1\left)\alpha h\right)& 1+n\alpha h\end{array}\right)$$

Then the linear system $\left(S_{h\alpha }\right)$ has a unique solution:

$${\mathrm{u}}_{\alpha ,i}=(b+\frac{g{x}_0-q}{\alpha })+(i-1)(g{x}_0-q)h-\frac{g}{\alpha }h-g\frac{i(i-1)}{2}{h}^2,i=1,\dots ,n+1.$$

As a consequence, the continuous solution $u_{\alpha }(x,y)$ of the system $\left(S_{\alpha }\right)$, can be approximated by the discrete function $u_{h\alpha }(x,y)$ in $\overline{\Omega }$ given by

$$u_{h\alpha }(x,y)=(g{x}_0-q-hgi)x+b+\frac{g{x}_0-q}{\alpha }-\frac{gh}{\alpha }+\frac{i^2-i}{2}g{h}^2,$$

(26)

for $x\in [x_i,x_{i+1}],y\in [0,y_0]$, $i=1,\dots ,n.$

Notice that ${\mathrm{u}}_{\alpha ,i}\to {\mathrm{u}}_i$, when $\alpha \to \infty$ for each $i=1,2,\dots ,n+1$. In addition, $u_{h\alpha }(x,y)\to u_h(x,y)$ when $\alpha \to \infty$ for every $(x,y)\in \overline{\Omega }$.

Lemma 2. 

If $u_{\alpha }$ is the solution of $\left(S_{\alpha }\right)$ and $u_{h\alpha }$ is the function given in (26) for each h, it follows that:

$${\parallel u_{\alpha }-u_{h\alpha }\parallel }_H\le C_{1\alpha }h,\mathit{and}{\parallel \frac{\partial u_{\alpha }}{\partial x}-\frac{\partial u_{h\alpha }}{\partial x}\parallel }_H\le \widetilde{C_1}h,$$

with $C_{1\alpha }=\left|g\right|x_0\sqrt{x_0{y}_0\left({\displaystyle \frac{2}{15}}+{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{\alpha x_0}}+\frac{1}{{\alpha }^2{x}_0^2}\right)}$ and $\widetilde{C_1}=\left|g\right|\sqrt{\frac{1}{3}{x}_0{y}_0}.$

Proof. 

From the definition of $u_{\alpha }$ and $u_{h\alpha }$, we have

$$\begin{array}{cc}& {\displaystyle \left|\right|u_{\alpha }-u_{h\alpha }{\left|\right|}_H^2=y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left(u_{\alpha }(x,y)-u_{h\alpha }(x,y)\right)}^{2}dx}\hfill \\ & {\displaystyle =y_0{g}^2\sum _{i=1}^n\frac{h^5{i}^2}{4}-\frac{h^{5}i}{6}+\frac{h^5}{20}-\frac{h^4}{3\alpha }+\frac{h^{4}i}{\alpha }+\frac{h^3}{{\alpha }^2}}\hfill \\ & =y_0{g}^2\left[h^5{n}^3\left({\displaystyle \frac{1}{12}}+{\displaystyle \frac{1}{24n}}+{\displaystyle \frac{1}{120n^2}}\right)+{\displaystyle \frac{h^4{n}^2}{\alpha }}\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{6n}}\right)+{\displaystyle \frac{h^{3}n}{{\alpha }^2}}\right]\hfill \\ & \le x_0{y}_0{g}^2{h}^2\left({\displaystyle \frac{2}{15}}{x}_0^2+{\displaystyle \frac{2}{3}}{\displaystyle \frac{x_0}{\alpha }}+{\displaystyle \frac{1}{{\alpha }^2}}\right)=C_{1\alpha }^2{h}^2.\hfill \end{array}$$

In addition, the derivatives of the functions $u_{\alpha }$ and $u_{h\alpha }$ are given by:

$$\frac{\partial u_{\alpha }}{\partial x}(x,y)=\frac{\partial u}{\partial x}(x,y)=-gx+g{x}_0-q,\forall (x,y)\in \Omega$$

and

$$\frac{\partial u_{h\alpha }}{\partial x}(x,y)=\frac{\partial u_h}{\partial x}(x,y)=g{x}_0-q-hgi,x\in [x_i,x_{i+1}],y\in [0,y_0].$$

Then, the bound for ${\parallel \frac{\partial u_{\alpha }}{\partial x}-\frac{\partial u_{h\alpha }}{\partial x}\parallel }_H$ coincides with the bound for ${\parallel \frac{\partial u}{\partial x}-\frac{\partial u_h}{\partial x}\parallel }_H$, obtained in Lemma 21. □

Remark 1. 

Notice that $C_{1\alpha }\to C_1$ when $\alpha \to \infty$, where $C_1$ is given in Lemma 21.

3. Distributed Optimization Problem with Variable g

In this section we are going to obtain discrete optimal solutions to the continuous optimization problem $\left(P_1\right)$ and $\left(P_{1\alpha }\right)$ in the rectangular domain $\Omega$ for the case when the optimization variable is g.

3.1. Discrete Problem $\left(P_{1h}\right)$ Associated to $\left(P_1\right)$

Taking into account that $b,g,q$ and the desired state $z_d$ are constants in the expression (6), according to [6], the continuous quadratic functional cost for the problem $\left(P_1\right)$ is explicitly given by:

$$J_1\left(g\right)=\frac{1}{2}{x}_0^3{y}_0{q}^2\left\{g^2\frac{x_0^2}{q^2}\left(\frac{2}{15}+\frac{M_1}{x_0^4}\right)+g\frac{x_0}{q}\left(-\frac{5}{12}+\frac{2}{3}\frac{(b-z_d)}{q{x}_0}\right)+\frac{1}{3}-\frac{(b-z_d)}{q{x}_0}+\frac{{(b-z_d)}^2}{q^2{x}_0^2}\right\}.$$

(27)

Since the function $J_1\left(g\right)$ is a polynomial function of the variable g, the analytic expresssion for $J_1^{\prime }\left(g\right)$ is given by

$$J_1^{\prime }\left(g\right)=\frac{1}{2}{x}_0^3{y}_0{q}^2\left\{\frac{2x_0^2}{q^2}g\left(\frac{2}{15}+\frac{M_1}{x_0^4}\right)+\frac{x_0}{q}\left(-\frac{5}{12}+\frac{2}{3}\frac{(b-z_d)}{q{x}_0}\right)\right\}$$

(28)

Then, the continuous solution to the distributed optimization problem $\left(P_1\right)$ is $g_{op}$ defined by:

$$g_{op}={\displaystyle \frac{q}{3x_0}}{\displaystyle \frac{\left(\frac{5}{8}-\frac{(b-z_d)}{q{x}_0}\right)}{\left(\frac{2}{15}+\frac{M_1}{x_0^4}\right)}}$$

(29)

and the continuous optimization state is

$$u_{g_{op}}(x,y)=-\frac{1}{2}{g}_{op}{x}^2+(g_{op}{x}_0-q)x+b.$$

(30)

We define the discrete distributed optimization problem $\left(P_{1h}\right)$ on the constant internal energy g as

$$\mathrm{find}{g}_{h_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{1h}\left(g_{h_{op}}\right)=\underset{g\in \mathbb{R}}{min}{J}_{1h}\left(g\right)$$

where the discrete cost function $J_{1h}$ is given by:

$$J_{1h}\left(g\right)=\frac{1}{2}{\parallel u_h-z_d\parallel }_H^2+\frac{1}{2}{M}_1{\parallel g\parallel }_H^2$$

with h defined in (14) and $u_h$ given in (22).

Taking into account that the variable g is constant it results that:

$$J_{1h}\left(g\right)=\frac{1}{2}{y}_0\left\{M_1{g}^2{x}_0+\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left(u_h(x,y)-z_d\right)}^{2}dx\right\}$$

(31)

and from algebraic work it follows that

$$\begin{array}{cc}\hfill J_{1h}\left(g\right)& ={\displaystyle \frac{1}{2}}{x}_0^3{y}_0{q}^2\{g^2\frac{x_0^2}{q^2}\left[{\displaystyle \frac{2}{15}}+\frac{M_1}{x_0^4}+{\displaystyle \frac{1}{180}}\frac{h^4}{x_0^4}+{\displaystyle \frac{1}{24}}\frac{h^3}{x_0^3}+{\displaystyle \frac{1}{36}}\frac{h^2}{x_0^2}-{\displaystyle \frac{5}{24}}\frac{h}{x_0}\right]\hfill \\ \hfill & +g\frac{x_0}{q}\left[-{\displaystyle \frac{5}{12}}+{\displaystyle \frac{2}{3}}\frac{(b-z_d)}{q{x}_0}+\frac{h}{x_0}\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{12}}\frac{h}{x_0}\right)-\frac{h}{x_0}\frac{(b-z_d)}{q{x}_0}\left({\displaystyle \frac{1}{2}}+\frac{1}{6}\frac{h}{x_0}\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{3}}-\frac{(b-z_d)}{q{x}_0}+\frac{{(b-z_d)}^2}{q^2{x}_0^2}\}.\hfill \end{array}$$

(32)

Lemma 3. 

The following estimate holds

$$|J_1\left(g\right)-J_{1h}\left(g\right)|\approx C_{2}h$$

(33)

where $C_2={\displaystyle \frac{1}{2}}{x}_0^3{y}_{0}gq|-{\displaystyle \frac{5}{24}}\frac{g{x}_0}{q}+\frac{1}{3}-\frac{1}{2}\frac{(b-z_d)}{q{x}_0}|$ does not depend on $h.$

Proof. 

From (27) and (32) we get

$$\begin{array}{cc}\hfill J_{1h}\left(g\right)-J_1\left(g\right)& ={\displaystyle \frac{1}{2}}{x}_0^3{y}_0{q}^2\{g^2\frac{x_0^2}{q^2}\left[{\displaystyle \frac{1}{180}}\frac{h^4}{x_0^4}+{\displaystyle \frac{1}{24}}\frac{h^3}{x_0^3}+{\displaystyle \frac{1}{36}}\frac{h^2}{x_0^2}-{\displaystyle \frac{5}{24}}\frac{h}{x_0}\right]\hfill \\ & +g\frac{h}{q}\left[{\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{12}}\frac{h}{x_0}-\frac{(b-z_d)}{q{x}_0}\left({\displaystyle \frac{1}{2}}+\frac{1}{6}\frac{h}{x_0}\right)\right]\}\hfill \\ & \approx {\displaystyle \frac{1}{2}}{x}_0^3{y}_0{q}^2\left[-{\displaystyle \frac{5}{24}}\frac{g^2{x}_0}{q^2}+\frac{g}{q}\left(\frac{1}{3}-\frac{1}{2}\frac{(b-z_d)}{q{x}_0}\right)\right]h\hfill \\ & ={\displaystyle \frac{1}{2}}{x}_0^3{y}_{0}gq\left[-{\displaystyle \frac{5}{24}}\frac{g{x}_0}{q}+\frac{1}{3}-\frac{1}{2}\frac{(b-z_d)}{q{x}_0}\right]h.\hfill \end{array}$$

Therefore, we obtain (33).

□

Since the function $J_{1h}\left(g\right)$ is a polynomial function of the variable g, it is easy to obtain the following analytic expression for the function $J_{1h}^{\prime }$:

$$\begin{array}{cc}\hfill J_{1h}^{\prime }\left(g\right)& ={\displaystyle \frac{1}{2}}{x}_0^3{y}_0{q}^2\{2g\frac{x_0^2}{q^2}\left[{\displaystyle \frac{2}{15}}+\frac{M_1}{x_0^4}+{\displaystyle \frac{1}{180}}\frac{h^4}{x_0^4}+{\displaystyle \frac{1}{24}}\frac{h^3}{x_0^3}+{\displaystyle \frac{1}{36}}\frac{h^2}{x_0^2}-{\displaystyle \frac{5}{24}}\frac{h}{x_0}\right]\hfill \\ \hfill & +\frac{x_0}{q}\left[-{\displaystyle \frac{5}{12}}+{\displaystyle \frac{2}{3}}\frac{(b-z_d)}{q{x}_0}+\frac{h}{x_0}\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{12}}\frac{h}{x_0}\right)-\frac{h}{x_0}\frac{(b-z_d)}{q{x}_0}\left({\displaystyle \frac{1}{2}}+\frac{1}{6}\frac{h}{x_0}\right)\right]\}.\hfill \end{array}$$

(34)

From the optimality condition we obtain the following result:

Lemma 4. 

a) 

The explicit expression for the optimal variable $g_{h_{op}}$ is given by:

$$g_{h_{op}}=\frac{q}{3x_0}\frac{A_1+\frac{h}{x_0}{A}_2+\frac{h^2}{x_0^2}{A}_3}{A_4+A_5\left(h\right)},$$

(35)

where

$$\begin{array}{cc}{A}_1=\frac{5}{8}-\frac{b-z_d}{q{x}_0},\hfill & A_4=\frac{2}{15}+\frac{M_1}{x_0^4},\hfill \\ A_2=\frac{3(b-z_d)}{4q{x}_0}-\frac{1}{2},\hfill & A_5\left(h\right)=\frac{h}{12x_0}\left(\frac{h^3}{15x_0^3}+\frac{h^2}{2x_0^2}+\frac{h}{3x_0}-\frac{5}{2}\right).\hfill \\ A_3=\frac{b-z_d}{4q{x}_0}-\frac{1}{8},\hfill \end{array}$$

(36)

b) 

In addition, the following error estimates hold:

$$|g_{op}-g_{h_{op}}|\approx C_{3}h,$$

(37)

$$\left|J_1\left(g_{op}\right)-J_{1h}\left(g_{h_{op}}\right)\right|\approx C_{4}h,$$

(38)

where $C_3$ and $C_4$ do not depend on h.

Proof. 

a)

It follows immediately from the expression of $J_{1h}^{\prime }\left(g\right)$ given by (34).

b)

Rewriting $g_{op}$ given by (29) as: $g_{op}=\frac{q}{3x_0}\frac{A_1}{A_4},$ it follows that:

$$g_{h_{op}}-g_{op}=\frac{q}{3x_0}\frac{-A_1{A}_5\left(h\right)+\frac{h}{x_0}{A}_2{A}_4+\frac{h^2}{x_0^2}{A}_3{A}_4}{A_4^2+A_4{A}_5\left(h\right)}\approx \frac{q}{3x_0^2}\frac{A_2{A}_4+\frac{5}{24}{A}_1}{A_4^2}h+o\left(h^2\right),$$

and we obtain (37) with $C_3=\left|C_3^{*}\right|$ where

$$C_3^{*}=\frac{q}{3x_0^2}\frac{A_2{A}_4+\frac{5}{24}{A}_1}{A_4^2}.$$

(39)

From the expression for $J_1\left(g\right)$ with $g=g_{op}$ and $J_{1h}\left(g\right)$ with $g=g_{h_{op}}$, it results that:

$$\begin{array}{cc}& J_1\left(g_{op}\right)-J_{1h}\left(g_{h_{op}}\right)=\frac{1}{2}{x}_0^3{y}_0{q}^2[\frac{x_0^2}{q^2}{A}_4(g_{op}^2-g_{h_{op}}^2)-\frac{2}{3}\frac{x_0}{q}{A}_1(g_{op}-g_{h_{op}})\hfill \\ & -g_{h_{op}}^2\frac{x_0^2}{q^2}{A}_5\left(h\right)+\frac{2}{3}{g}_{h_{op}}\frac{h}{q}\left(A_2+\frac{h}{x_0}{A}_3\right)].\hfill \end{array}$$

By using (37) we get (38) where

$$\begin{array}{cc}{C}_4\hfill & =\frac{1}{2}{x}_0^2{y}_0{q}^2\left|-2A_4{C}_3^{*}{g}_{op}\frac{x_0^3}{q^2}+\frac{5}{24}{g}_{op}^2\frac{x_0^2}{q^2}+\frac{2}{3}{A}_1{C}_3^{*}\frac{x_0^2}{q}+\frac{2}{3}{A}_2{g}_{op}\frac{x_0}{q}\right|\hfill \\ \\ & =\frac{1}{2}{x}_0^2{y}_0{q}^2\left|\frac{A_1(5A_1+48A_2{A}_4)}{216A_4^2}\right|.\hfill \end{array}$$

□

Lemma 5. 

Let us consider $u_{g_{op}}$ the solution of (1)-(2) for $g=g_{op}$ and $u_{h{g}_{h_{op}}}$ the discrete solution defined as in (22) for $h>0$, where $g=g_{h_{op}}$ is the optimal variable of $\left(P_{1h}\right)$ given by (35). We have that:

$$a){\parallel u_{g_{op}}-u_{h{g}_{h_{op}}}\parallel }_H\approx C_{5}h,b){\parallel \frac{\partial u_{g_{op}}}{\partial x}-\frac{\partial u_{h{g}_{h_{op}}}}{\partial x}\parallel }_H\approx C_{6}h$$

where $C_5$ and $C_6$ are positive constants independent of the parameter $h.$

Proof:

a)

From definition of norm in H, we obtain

$$\begin{array}{cc}& {\displaystyle \left|\right|u_{g_{op}}-u_{h{g}_{h_{op}}}{\left|\right|}_H^2=y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left(u_{g_{op}}(x,y)-u_{h{g}_{h_{op}}}(x,y)\right)}^{2}dx}\hfill \\ & {\displaystyle =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left[-\frac{1}{2}{g}_{op}{x}^2+x_{0}x(g_{op}-g_{h_{op}})+h{g}_{h_{op}}(ix-h\frac{i(i-1)}{2})\right]}^{2}dx}\hfill \\ & =\frac{1}{120}{x}_0^3{y}_0{g}_{op}^2\left[10-25\frac{x_0{C}_3^{*}}{g_{op}}+16{\left(\frac{x_0{C}_3^{*}}{g_{op}}\right)}^2\right]h^2+o\left(h^3\right).\hfill \end{array}$$

where $C_3^{*}$ is given by (39). Therefore it follows that

$${\parallel u_{g_{op}}-u_{h{g}_{h_{op}}}\parallel }_H\approx C_{5}h,$$

with

$$C_5=\left|g_{op}\right|\sqrt{\frac{1}{120}{x}_0^3{y}_0\left[10-25\frac{x_0{C}_3^{*}}{g_{op}}+16{\left(\frac{x_0{C}_3^{*}}{g_{op}}\right)}^2\right]}.$$

b)

We have that

$\begin{array}{cc}& {∥\frac{\partial u_{g_{op}}}{\partial x}-\frac{\partial u_{h{g}_{h_{op}}}}{\partial x}∥}_H^2\hfill \\ & {\displaystyle =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{\left(-g_{op}x+h{g}_{h_{op}}i+x_0(g_{op}-g_{h_{op}})\right)}^{2}dx}\hfill \\ & =\frac{x_0{y}_0}{6}\left[2g_{op}^2{x}_0^2+g_{op}{g}_{h_{op}}(h^2+3h{x}_0-4x_0^2)+g_{h_{op}}^2(h^2-3h{x}_0+2x_0^2)\right]\hfill \\ & =\frac{x_0{y}_0}{6}\left[2x_0^2{(g_{op}-g_{h_{op}})}^2+3h{x}_0{g}_{h_{op}}(g_{op}-g_{h_{op}})+h^2{g}_{h_{op}}(g_{op}+g_{h_{op}})\right].\hfill \end{array}$

Taking into account Lemma 4, we get

$${∥\frac{\partial u_{g_{op}}}{\partial x}-\frac{\partial u_{h{g}_{o{p}_h}}}{\partial x}∥}_H^2=\frac{x_0{y}_0}{6}{g}_{op}^2\left[2+3\frac{x_0{C}_3^{*}}{g_{op}}+2{\left(\frac{x_0{C}_3^{*}}{g_{op}}\right)}^2\right]h^2+o\left(h^3\right).$$

Then

$${∥\frac{\partial u_{g_{op}}}{\partial x}-\frac{\partial u_{h{g}_{o{p}_h}}}{\partial x}∥}_H\approx C_{6}h,$$

with

$$C_6=\left|g_{op}\right|\sqrt{\frac{x_0{y}_0}{6}\left[2-3\frac{x_0{C}_3^{*}}{g_{op}}+2{\left(\frac{x_0{C}_3^{*}}{g_{op}}\right)}^2\right]},$$

where $C_3^{*}$ is given by (39).

3.2. Discrete Problem $\left(P_{1h\alpha }\right)$ Associated to $\left(P_{1\alpha }\right)$

From [6], we know that the continuous quadratic functional cost in (6) for the optimization problem $\left(P_{1\alpha }\right)$ is explicitly given by:

$$\begin{array}{c}{J}_{1\alpha }\left(g\right)=J_1\left(g\right)\hfill \\ +\frac{x_0^2{y}_0{q}^2}{2\alpha }\left\{\frac{g^2{x}_0^2}{q^2}\left(\frac{2}{3}+\frac{1}{\alpha x_0}\right)+\frac{g{x}_0}{q}\left(-\frac{5}{3}-\frac{2}{\alpha x_0}+\frac{2(b-z_d)}{q{x}_0}\right)+1+\frac{1}{\alpha x_0}-\frac{2(b-z_d)}{q{x}_0}\right)\},\hfill \end{array}$$

(40)

where $J_1$ is defined by (27). Moreover, the continuous optimal distributed variable denoted by $g_{{\alpha }_{op}}$ is

$$g_{{\alpha }_{op}}={\displaystyle \frac{q}{3x_0}}{\displaystyle \frac{\left(\frac{5}{8}-\frac{(b-z_d)}{q{x}_0}+\frac{5}{2\alpha x_0}+\frac{3}{{\alpha }^2{x}_0^2}-\frac{3(b-z_d)}{\alpha q{x}_0^2}\right)}{\left(\frac{2}{15}+\frac{M_1}{x_0^4}+\frac{2}{3\alpha x_0}+\frac{1}{{\alpha }^2{x}_0^2}\right)}}.$$

(41)

The continuous associated state is established by:

$$u_{\alpha g_{{\alpha }_{op}}}(x,y)=-\frac{1}{2}{g}_{{\alpha }_{op}}{x}^2+(g_{{\alpha }_{op}}{x}_0-q)x+\frac{1}{\alpha }(g_{{\alpha }_{op}}{x}_0-q)+b.$$

(42)

Defining the discrete cost function as:

$$J_{1h\alpha }\left(g\right)=\frac{1}{2}{\parallel u_{h\alpha }-z_d\parallel }_H^2+\frac{1}{2}{M}_1{\parallel g\parallel }_H^2$$

(43)

where the function $u_{h\alpha }$ is given in (26) we set the following discrete optimization problem $\left(P_{1h\alpha }\right)$ on the constant internal energy g as

$$\mathrm{find}{g}_{h{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{1h\alpha }\left(g_{h{\alpha }_{op}}\right)=\underset{g\in \mathbb{R}}{min}{J}_{1\alpha h}\left(g\right).$$

The discrete cost function $J_{1h\alpha }$ is explicitly given by

$$\begin{array}{cc}\hfill & J_{1h\alpha }\left(g\right)=J_{1\alpha }\left(g\right)\hfill \\ \hfill & +\frac{1}{2}{x}_0^3{y}_{0}gqh\left\{\frac{g{x}_0}{q}\right[-\frac{5}{24}+\frac{1}{36}\frac{h}{x_0}+\frac{1}{24}\frac{h^2}{x_0^2}+\frac{1}{180}\frac{h^3}{x_0^3}+\frac{1}{\alpha x_0}\left(-\frac{7}{6}+\frac{1}{3}\frac{h}{x_0}+\frac{1}{6}\frac{h^2}{x_0^2}\right)\hfill \\ \hfill & +\frac{1}{{\alpha }^2{x}_0^2}\left(-2+\frac{h}{x_0}\right)]+\frac{1}{3}+\frac{1}{12}\frac{h}{x_0}+\frac{1}{\alpha x_0}\left(\frac{3}{2}+\frac{1}{6}\frac{h}{x_0}\right)+\frac{2}{{\alpha }^2{x}_0^2}\hfill \\ \hfill & +\frac{(b-z_d)}{q{x}_0}\left(-\frac{1}{2}-\frac{1}{6}\frac{h}{x_0}-\frac{2}{\alpha x_0}\right)\},\hfill \end{array}$$

(44)

where $J_{1\alpha }$ is given by (40).

Lemma 6. 

For $g\in H,h>0$ and $\alpha >0$ the following estimate holds

$$|J_{1h\alpha }\left(g\right)-J_{1\alpha }\left(g\right)|\approx C_{2\alpha }h$$

with

$$C_{2\alpha }=\frac{1}{2}{x}_0^3{y}_0\left|g\right||\frac{g{x}_0}{q}\left(-\frac{5}{24}-\frac{7}{6}\frac{1}{\alpha x_0}-\frac{2}{{\alpha }^2{x}_0^2}\right)+\frac{1}{3}+\frac{3}{2}\frac{1}{\alpha x_0}+\frac{2}{{\alpha }^2{x}_0^2}+\frac{(b-z_d)}{q{x}_0}\left(-\frac{1}{2}-\frac{2}{\alpha x_0}\right)|,$$

a constant independent of $h.$

Proof. 

It follows immediately from expression (44). □

Remark 2. 

$C_{2\alpha }\to C_2$ when $\alpha \to \infty$, where $C_2$ is given in Lemma 3.

Lemma 7. 

a) 

The explicit expression for the optimal variable $g_{h{\alpha }_{op}}$ is given by:

$$g_{h{\alpha }_{op}}=\frac{q}{3x_0}\frac{A_{1\alpha }+\frac{h}{x_0}{A}_{2\alpha }+\frac{h^2}{x_0^2}{A}_{3\alpha }}{A_{4\alpha }+A_{5\alpha }\left(h\right)},$$

(45)

where

$$\begin{array}{cc}\hfill A_{1\alpha }& =\frac{5}{8}-\frac{b-z_d}{q{x}_0}+\frac{1}{\alpha x_0}\left(\frac{5}{2}+\frac{3}{\alpha x_0}-\frac{3(b-z_d)}{q{x}_0}\right),\hfill \\ \hfill A_{2\alpha }& =\frac{3(b-z_d)}{4q{x}_0}-\frac{1}{2}+\frac{1}{\alpha x_0}\left(-\frac{9}{4}-\frac{3}{\alpha x_0}+\frac{3(b-z_d)}{q{x}_0}\right),\hfill \\ \hfill A_{3\alpha }& =\frac{b-z_d}{4q{x}_0}-\frac{1}{8}-\frac{1}{4\alpha x_0},\hfill \\ \hfill A_{4\alpha }& =\frac{2}{15}+\frac{M_1}{x_0^4}+\frac{1}{\alpha x_0}\left(\frac{2}{3}+\frac{1}{\alpha x_0}\right),\hfill \\ \hfill A_{5\alpha }\left(h\right)& =\frac{h}{12x_0}\left(\frac{h^3}{15x_0^3}+\frac{h^2}{2x_0^2}+\frac{h}{3x_0}-\frac{5}{2}\right)\hfill \\ \hfill & +\frac{h}{\alpha x_0^2}\left(-\frac{7}{6}+\frac{h}{3x_0}+\frac{h^2}{6x_0^2}+\frac{1}{\alpha x_0}\left(-2+\frac{h}{x_0}\right)\right).\hfill \end{array}$$

(46)

b) 

In addition, the following error estimates hold:

$$|g_{{\alpha }_{op}}-g_{h{\alpha }_{op}}|\approx C_{3\alpha }h,$$

(47)

$$\left|J_{1\alpha }\left(g_{{\alpha }_{op}}\right)-J_{1h\alpha }\left(g_{h{\alpha }_{op}}\right)\right|\approx C_{4\alpha }h,$$

(48)

where $C_{3\alpha }$ and $C_{4\alpha }$ do not depend on h.

Proof. 

a)

It follows immediately from the fact that

$$J_{1h\alpha }\left(g\right)=\frac{1}{2}{x}_0^3{y}_0{q}^2\left[2g\frac{x_0^2}{q^2}\left(A_{4\alpha }+A_{5\alpha }\left(h\right)\right)-\frac{2}{3}\frac{x_0}{q}\left(A_{1\alpha }+\frac{h}{x_0}{A}_{2\alpha }+\frac{h^2}{x_0^2}{A}_{3\alpha }\right)\right].$$

b)

Notice that $g_{{\alpha }_{op}}$ given by (41) can be rewriten as $g_{{\alpha }_{op}}=\frac{q}{3x_0}\frac{A_{1\alpha }}{A_{4\alpha }}.$ Then,

$$\begin{array}{cc}{g}_{h{\alpha }_{op}}-g_{{\alpha }_{op}}\hfill & =\frac{q}{3x_0}\frac{-A_{1\alpha }{A}_{5\alpha }\left(h\right)+\frac{h}{x_0}{A}_{2\alpha }{A}_{4\alpha }+\frac{h^2}{x_0^2}{A}_{3\alpha }{A}_{4\alpha }}{A_{4\alpha }^2+A_{4\alpha }{A}_{5\alpha }\left(h\right)}\hfill \\ & \approx \frac{q}{3x_0^2}\frac{A_{2\alpha }{A}_{4\alpha }+\frac{5}{24}{A}_{1\alpha }}{A_{4\alpha }^2}h+o\left(h^2\right).\hfill \end{array}$$

and we obtain (47) with $C_{3\alpha }=\left|C_{3\alpha }^{*}\right|$ where

$$C_{3\alpha }^{*}=\frac{q}{3x_0^2}\frac{A_{2\alpha }{A}_{4\alpha }+\left(\frac{5}{24}+\frac{7}{6\alpha x_0}\right)A_{1\alpha }}{A_{4\alpha }^2}.$$

(49)

Following Lemma 4 it is obtained formula (48) with

$$\begin{array}{cc}{C}_{4\alpha }\hfill & =\frac{1}{2}{x}_0^2{y}_0{q}^2|-2A_{4\alpha }{C}_{3\alpha }^{*}{g}_{{\alpha }_{op}}\frac{x_0^3}{q^2}+\frac{5}{24}{g}_{{\alpha }_{op}}^2\frac{x_0^2}{q^2}\hfill \\ & +\frac{7}{6}{g}_{{\alpha }_{op}}^2\frac{x_0}{\alpha q^2}+\frac{2}{3}{A}_{1\alpha }{C}_{3\alpha }^{*}\frac{x_0^2}{q}+\frac{2}{3}{A}_{2\alpha }{g}_{{\alpha }_{op}}\frac{x_0}{q}|.\hfill \end{array}$$

□

Remark 3. 

When $\alpha \to \infty$, we have that $A_{i\alpha }\to A_i$, where $A_i$ and $A_{i\alpha }$ are given by (36) and (46), respectively for $i=1,2,\dots ,5$. As an immediate consequence it follows that, $C_{3\alpha }\to C_3$ and $C_{4\alpha }\to C_4$ when $\alpha \to \infty$, where $C_3$ and $C_4$ are defined in Lemma 4.

Lemma 8. 

Let us consider $u_{\alpha g_{{\alpha }_{op}}}$ the function given by (14) for $g=g_{{\alpha }_{op}}$ where $g_{{\alpha }_{op}}$ is the optimal variable of problem $\left(P_{1\alpha }\right)$ given by (41) and $u_{h\alpha g_{h{\alpha }_{op}}}$ the function defined by (26) for $h>0$ where $g=g_{h{\alpha }_{op}}$ is the optimal variable of $\left(P_{1h\alpha }\right)$ given by (35). We have that:

$$a){\parallel u_{\alpha g_{{\alpha }_{op}}}-u_{h\alpha g_{h{\alpha }_{op}}}\parallel }_H\approx C_{5\alpha }h,b){∥\frac{\partial u_{\alpha g_{{\alpha }_{op}}}}{\partial x}-\frac{\partial u_{h\alpha g_{h{\alpha }_{op}}}}{\partial x}∥}_H\approx C_{6\alpha }h,$$

where $C_{5\alpha }$ and $C_{6\alpha }$ are positive constants independent of the parameter $h.$

Proof. 

Working algebraically we can obtain that

$$\begin{array}{c}{C}_{5\alpha }=\left|g_{{\alpha }_{op}}\right|\left\{\frac{1}{120}{x}_0^3{y}_0\right[10-25\frac{x_0{C}_{3\alpha }^{*}}{g_{op}}+16{\left(\frac{x_0{C}_{3\alpha }^{*}}{g_{op}}\right)}^2\hfill \\ +\frac{1}{\alpha x_0}\left(60+\frac{120}{\alpha x_0}-240\frac{C_{3\alpha }^{*}}{\alpha g_{{\alpha }_{op}}}+120\frac{C_{3\alpha }^{*}{x}_0}{\alpha g_{{\alpha }_{op}}^2}-140\frac{C_{3\alpha }^{*}{x}_0}{g_{{\alpha }_{op}}}+80\frac{C_{3\alpha }^{*}{x}_0^2}{g_{{\alpha }_{op}}^2}\right)]{\}}^{1/2},\hfill \end{array}$$

and

$$C_{6\alpha }=\left|g_{{\alpha }_{op}}\right|\sqrt{\frac{x_0{y}_0}{6}\left[2-3\frac{x_0{C}_{3\alpha }^{*}}{g_{{\alpha }_{op}}}+2{\left(\frac{x_0{C}_{3\alpha }^{*}}{g_{{\alpha }_{op}}}\right)}^2\right]},$$

where $C_{3\alpha }^{*}$ is given by (49). □

Remark 4. 

$C_{5\alpha }\to C_5$ and $C_{6\alpha }\to C_6$ when $\alpha \to \infty$ where $C_5$ and $C_6$ are given in Lemma 5.

Remark 5. 

In [23] the double convergence when $(h,\alpha )\to (0,+\infty )$ of the optimal control problem $\left(P_{1h\alpha }\right)$ has been studied, obtaining a commutative diagram that relates the continuous and discrete optimal control problems $\left(P_1\right)$, $\left(P_{1\alpha }\right)$, $\left(P_{1h}\right)$ and $\left(P_{1h\alpha }\right)$ as in the following scheme:

4. Boundary Optimization Problem with Variable q

4.1. Discrete Problem $\left(P_{2h}\right)$ Associated to $\left(P_2\right)$

Under the same considerations given in Section 3.1 and taking into account formula (9), for a given $q\in Q$, we obtain the following quadratic cost function:

$$\begin{array}{cc}{J}_2\left(q\right)\hfill & =\frac{x_0{y}_0}{2}\{q^2{x}_0^2\left(\frac{1}{3}+\frac{M_2}{x_0^3}\right)+q{x}_0\left(-\frac{5}{12}g{x}_0^2-(b-z_d)\right)\hfill \\ \hfill & +\frac{2}{15}{g}^2{x}_0^4+{(b-z_d)}^2+\frac{2}{3}g{x}_0^2(b-z_d)\}.\hfill \end{array}$$

(50)

Then, the boundary optimal control of the problem $\left(P_2\right)$, called $q_{op}$, and the associated continuous optimal state are given by:

$$q_{op}={\displaystyle \frac{\frac{5}{12}g{x}_0^2+(b-z_d)}{2x_0\left(\frac{1}{3}+\frac{M_2}{x_0^3}\right)}},u_{q_{op}}(x,y)=-\frac{1}{2}g{x}^2+(g{x}_0-q_{op})x+b.$$

(51)

Associated to $\left(P_2\right)$, we define the approximate discrete distributed optimal control problem $\left(P_{2h}\right)$ on the constant heat flux q as

$$\mathrm{find}{q}_{h_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{2h}\left(q_{h_{op}}\right)=\underset{q\in \mathbb{R}}{min}{J}_{2h}\left(q\right)$$

where the discrete cost function $J_{2h}$ is defined by

$$J_{2h}\left(q\right)=\frac{1}{2}{\parallel u_h-z_d\parallel }_H^2+\frac{1}{2}{M}_2{\parallel q\parallel }_Q^2,$$

where $u_h$ is given in (22), h is the spatial step, and $z_d$ (the desired state) and the control q are constant. From the definition of norm over Q, it results that:

$$J_{2h}\left(q\right)=\frac{1}{2}{y}_0\left\{M_2{q}^2{x}_0+\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{[u_h(x,y)-z_d]}^{2}dx\right\}$$

and working algebraically, we get

$$\begin{array}{cc}\hfill & J_{2h}\left(q\right)=J_2\left(q\right)+\frac{x_0{y}_0}{2}\{hg{x}_0\left[\frac{1}{3}q{x}_0-\frac{5}{24}g{x}_0^2-\frac{1}{2}(b-z_d)\right]\hfill \\ \hfill & +h^{2}g\left[\frac{1}{36}{x}_0^{2}g-\frac{1}{6}(b-z_d)+\frac{1}{12}q{x}_0\right]+\frac{1}{24}{h}^3{g}^2{x}_0+\frac{1}{180}{h}^4{g}^2\}.\hfill \end{array}$$

(52)

Lemma 9. 

Given $q\in Q$ and $h>0$, we have that

$$|J_{2h}\left(q\right)-J_2\left(q\right)|\approx C_{7}h,$$

(53)

where $C_7=\frac{1}{2}{x}_0^2{y}_0\left|g\right|\left|\frac{1}{3}q{x}_0-\frac{5}{24}g{x}_0^2-\frac{1}{2}(b-z_d)\right|$ is a constant independent of $h.$

proof: It follows immediately from the expression (52) for $J_{2h}$.

Lemma 10. 

Let us consider $h>0.$

a)

The explicit expression for the optimal variable $q_{h_{op}}$ is given by:

$$q_{h_{op}}=q_{op}-B_1\frac{gh}{6}-B_2\frac{g{h}^2}{24x_0},$$

(54)

with $B_1=B_2={\left(\frac{1}{3}+\frac{M_2}{x_0^3}\right)}^{-1}.$

b)

The following error estimates hold:

$$|q_{h_{op}}-q_{op}|\approx C_{8}h,$$

(55)

$$\left|J_{2h}\left(q_{h_{op}}\right)-J_2\left(q_{op}\right)\right|\approx C_{9}h,$$

(56)

where $C_8$ and $C_9$ are constants independent of h.

Proof. 

a)

From the expression (52) for $J_{2h}$ we have that

$$J_{2h}^{\prime }\left(q\right)=\frac{1}{2}{x}_0^2{y}_0\left\{2x_{0}q\left(\frac{1}{3}+\frac{M_2}{x_0^3}\right)-\frac{5}{12}g{x}_0^2-(b-z_d)+\frac{1}{3}g{x}_{0}h+\frac{1}{12}g{h}^2\right\}.$$

Therefore

$$q_{h_{op}}={\displaystyle \frac{\frac{5}{12}g{x}_0^2+(b-z_d)-\frac{1}{12}g{h}^2-\frac{1}{3}g{x}_{0}h}{2x_0\left(\frac{1}{3}+\frac{M_2}{x_0^3}\right)}}.$$

Taking into account that $q_{op}$ is given by (51), we obtain formula (54).

b)

On one hand, expression (55) is a direct consequence of expression (54) with $C_8=\left|\frac{g{B}_1}{6}\right|.$

On the other hand, taking into account formula (52) for $J_{2h}$ it follows that

$$\begin{array}{c}{J}_{2h}\left(q_{h_{op}}\right)-J_2\left(q_{op}\right)=J_2\left(q_{h_{op}}\right)-J_2\left(q_{op}\right)\hfill \\ +\frac{x_0^2{y}_0}{2}hg\left[\frac{1}{3}{q}_{h_{op}}{x}_0-\frac{5}{24}g{x}_0^2-\frac{1}{2}(b-z_d)\right]+o\left(h^2\right).\hfill \end{array}$$

(57)

From the definition of $J_2$ given by (50) and the explicit expression (54) for $q_{h_{op}}$ we get

$$\begin{array}{cc}\hfill & J_2\left(q_{h_{op}}\right)-J_2\left(q_{op}\right)=\frac{x_0^2{y}_0}{2}\{\frac{x_0}{B_1}(q_{h_{op}}^2-q_{op}^2)\hfill \\ \hfill & +(q_{h_{op}}-q_{op})\left(-\frac{5}{12}g{x}_0^2-(b-z_d)\right)\}\hfill \\ \hfill & =\frac{x_0^2{y}_0}{2}(q_{h_{op}}-q_{op})\left(\frac{x_0}{B_1}(q_{h_{op}}+q_{op})-\frac{5}{12}g{x}_0^2-(b-z_d)\right)\hfill \\ \hfill & =\frac{x_0^2{y}_0}{12}gh{B}_1\left(-\frac{2x_0}{B_1}{q}_{op}+\frac{5}{12}g{x}_0^2+b-z_d\right)+o\left(h^2\right).\hfill \end{array}$$

(58)

Taking into account (57) and (58), we obtain the estimate (56) with

$$C_9={\displaystyle \frac{\left|g\right|x_0^2{y}_0}{2}}\left|\frac{B_1}{6}\left(-\frac{2x_0}{B_1}{q}_{op}+\frac{5}{12}g{x}_0^2+b-z_d\right)+\frac{1}{3}{q}_{op}{x}_0-\frac{5}{24}g{x}_0^2-\frac{1}{2}(b-z_d)\right|.$$

□

Lemma 11. 

Consider $u_{q_{op}}$ the solution of (1)-(2) for $q=q_{op}$ and $u_{h{q}_{h_{op}}}$ the discrete solution given by (22) for each $h>0$ where $q=q_{h_{op}}$ is the optimal variable of the problem $\left(P_{2h}\right)$ given by (54). Then we have that:

$$\mathrm{a}){\parallel u_{h{q}_{h_{op}}}-u_{q_{op}}\parallel }_H\approx C_{10}h,\mathrm{b}){∥\frac{\partial u_{h{q}_{h_{op}}}}{\partial x}-\frac{\partial u_{q_{op}}}{\partial x}∥}_H\approx C_{11}h,$$

(59)

where $C_{10}$ and $C_{11}$ do not depend on the parameter $h.$

Proof. 

a)

From the definition of u and $u_h$ given by (13) and (22), respectively, it follows for $i=1,\dots ,n$ that

$$\begin{array}{c}{u}_{h{q}_{h_{op}}}(x,y)-u_{q_{op}}(x,y)=\frac{1}{2}g{x}^2-(q_{h_{op}}-q_{op}+gih)x+g\frac{i^2-i}{2}{h}^2,\hfill \\ =\frac{1}{2}g{x}^2+hgx\left(\frac{B_1}{6}-i\right)+h^{2}g\left(\frac{B_1}{24}\frac{x}{x_0}+\frac{i(i-1)}{2}\right),\forall x\in [x_i,x_{i+1}].\hfill \end{array}$$

(60)

Then

$$\begin{array}{c}{\displaystyle {\parallel u_{h{q}_{h_{op}}}-u_{q_{op}}\parallel }_H^2=y_0\sum _{i=1}^n\underset{x_i}{\overset{x_{i+1}}{\int }}{\left(u_{h{q}_{h_{op}}}\left(x\right)-u_{q_{op}}\left(x\right)\right)}^{2}dx}\hfill \\ =x_0{y}_0{g}^2[\frac{1}{108}{(B_1-3)}^2{h}^2{x}_0^2+\frac{1}{216}{h}^3{x}_0{(B_1-3)}^2\hfill \\ +{\displaystyle \frac{1}{8640}}{h}^4(72-30B_1+5B_1^2)].\hfill \end{array}$$

Therefore, we get formula (59) a) with $C_{10}=\left|g(B_1-3)\right|x_0\sqrt{\frac{x_0{y}_0}{108}}$ .

b)

In the same manner, we get

$$\begin{array}{cc}& {\displaystyle {∥\frac{\partial u_{h{q}_{h_{op}}}}{\partial x}-\frac{\partial u_{q_{op}}}{\partial x}∥}_H^2=y_0\sum _{i=1}^n\underset{x_i}{\overset{x_{i+1}}{\int }}{\left(gx-(q_{o{p}_h}-q_{op})-hgi\right)}^{2}dx}\hfill \\ & =y_0{g}^2\left[\frac{1}{36}(12-6B_1+B_1^2)h^2{x}_0+\frac{1}{72}(B_1-3)B_1{h}^3+\frac{1}{576}{B}_1^2{\displaystyle \frac{h^4}{x_0}}\right].\hfill \end{array}$$

Then we get (59) b) with

$$C_{11}={\displaystyle \frac{\left|g\right|}{6}}\sqrt{x_0{y}_0(12-6B_1+B_1^2)}.$$

□

4.2. Discrete Problem $\left(P_{2h\alpha }\right)$ Associated to $\left(P_{2\alpha }\right)$

If we suppose that the desired state $z_d$ is constant in (9), the quadratic cost function $J_{2\alpha }$ for the optimal control problem $\left(P_{2\alpha }\right)$ is explicitly given by:

$$\begin{array}{c}{J}_{2\alpha }\left(q\right)={\displaystyle \frac{x_0{y}_0}{2}}[q^2{x}_0^2(D_{1\alpha }+\frac{M_2}{x_0^3})+q{x}_0\left(D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-z_d)\right)\hfill \\ +D_{4\alpha }{g}^2{x}_0^4+D_{5\alpha }{(b-z_d)}^2+D_{6\alpha }g{x}_0^2(b-z_d)],\hfill \end{array}$$

(61)

where

$$\begin{array}{cc}{D}_{1\alpha }=\frac{1}{3}+\frac{1}{\alpha x_0}+\frac{1}{{\alpha }^2{x}_0^2},\hfill & D_{2\alpha }=-\frac{5}{12}-\frac{5}{3\alpha x_0}-\frac{2}{{\alpha }^2{x}_0^2},\hfill \\ D_{3\alpha }=-1-\frac{2}{\alpha x_0},\hfill & D_{4\alpha }=\frac{2}{15}+\frac{2}{3\alpha x_0}+\frac{1}{{\alpha }^2{x}_0^2},\hfill \\ D_{5\alpha }=1,\hfill & D_{6\alpha }=\frac{2}{3}+\frac{2}{\alpha x_0}.\hfill \end{array}$$

(62)

Then, the continuous boundary optimization control, called $q_{{\alpha }_{op}}$, and the associated state are:

$$\begin{array}{c}{q}_{{\alpha }_{op}}=-\frac{D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-z_d)}{2x_0\left(D_{1\alpha }+\frac{M_2}{x_0^3}\right)},\hfill \\ \\ u_{\alpha q_{{\alpha }_{op}}}(x,y)=-\frac{1}{2}g{x}^2+(g{x}_0-q_{{\alpha }_{op}})x+\frac{1}{\alpha }(g{x}_0-q_{{\alpha }_{op}})+b.\hfill \end{array}$$

(63)

Remark 6. 

Notice that $J_{2\alpha }\left(q\right)\to J_2\left(q\right)$ for all $q\in Q$ and $q_{{\alpha }_{op}}\to q_{op}$ when $\alpha \to \infty$.

Defining the discrete cost function as:

$$J_{2h\alpha }\left(q\right)=\frac{1}{2}{\parallel u_{h\alpha }-z_d\parallel }_H^2+\frac{1}{2}{M}_2{\parallel q\parallel }_Q^2,$$

(64)

where $u_{h\alpha }$ is the solution of $\left(S_{h\alpha }\right)$ given in (26), we set the following discrete optimization problem $\left(P_{2h\alpha }\right)$ on the constant heat flux q as

$$\mathrm{find}{q}_{h{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{2h\alpha }\left(q_{h{\alpha }_{op}}\right)=\underset{q\in \mathbb{R}}{min}{J}_{2h\alpha }\left(g\right).$$

Working algebraically, the cost function $J_{2h\alpha }$ can be written explicitly as:

$$\begin{array}{cc}\hfill & J_{2h\alpha }\left(q\right)=J_{2\alpha }\left(q\right)\hfill \\ \hfill & +\frac{1}{2}{x}_0{y}_{0}gh\{\frac{q{x}_0^2}{3}-\frac{5g{x}_0^3}{24}-\frac{(b-z_d)x_0}{2}+\frac{3x_{0}q}{2\alpha }-\frac{7g{x}_0^2}{6\alpha }-\frac{2(b-z_d)}{\alpha }+\frac{2(q-g{x}_0)}{{\alpha }^2}\hfill \\ \hfill & +h\left(\frac{x_0^{2}g}{36}-\frac{(b-z_d)}{6}+\frac{q{x}_0}{12}+\frac{2g{x}_0+q}{6\alpha }+\frac{g}{{\alpha }^2}\right)\hfill \\ \hfill & +h^{2}g\left(\frac{x_0}{24}+\frac{1}{6\alpha }\right)+\frac{1}{180}g{h}^3\}.\hfill \end{array}$$

(65)

Lemma 12. 

For each $q\in Q$ and $h>0$, we have that:

$$|J_{2h\alpha }\left(q\right)-J_{2\alpha }\left(q\right)|\approx C_{7\alpha }h,$$

(66)

with

$$C_{7\alpha }=\frac{1}{2}{x}_0^2{y}_0\left|g\left(\frac{q{x}_0}{3}-\frac{5g{x}_0^2}{24}-\frac{(b-z_d)}{2}+\frac{3q}{2\alpha }-\frac{7g{x}_0}{6\alpha }-\frac{2(b-z_d)}{x_0\alpha }+\frac{2(q-g{x}_0)}{x_0{\alpha }^2}\right)\right|,$$

a constant independent of $h.$

Proof. 

It follows immediately from expression (65). □

Lemma 13. 

Let us consider $h>0.$

a)

The explicit expression for the optimal control $q_{h_{op}}$ is given by:

$$q_{h{\alpha }_{op}}=q_{{\alpha }_{op},}-B_{1\alpha }{\displaystyle \frac{gh}{6}}-B_{2\alpha }{\displaystyle \frac{g{h}^2}{24x_0}}$$

(67)

with

$$\begin{array}{c}{B}_{1\alpha }=\left(1+\frac{9}{2\alpha x_0}+\frac{6}{{\alpha }^2{x}_0^2}\right){\left(D_{1\alpha }+\frac{M_2}{x_0^3}\right)}^{-1},\\ B_{2\alpha }=\left(1-\frac{2}{\alpha x_0}\right){\left(D_{1\alpha }+\frac{M_2}{x_0^3}\right)}^{-1}.\end{array}$$

b)

The following error estimates hold:

$$|q_{h{\alpha }_{op}}-q_{{\alpha }_{op}}|\approx C_{8\alpha }h,$$

(68)

$$\left|J_{2h\alpha }\left(q_{h{\alpha }_{op}}\right)-J_{2\alpha }\left(q_{{\alpha }_{op}}\right)\right|\approx C_{9\alpha }h,$$

(69)

where $C_{8\alpha }$ and $C_{9\alpha }$ does not depend on h.

Proof. 

a)

From the derivative of the control function $J_{2h\alpha }$ given by

$$\begin{array}{c}{J}_{2h\alpha }^{\prime }\left(q\right)=J_{2\alpha }^{\prime }\left(q\right)+\frac{1}{2}{x}_0{y}_{0}gh\left\{\frac{x_0^2}{3}+\frac{3x_0}{2\alpha }+\frac{2}{{\alpha }^2}+h\left(\frac{x_0}{12}+\frac{1}{6\alpha }\right)\right\}\hfill \\ ={\displaystyle \frac{x_0{y}_0}{2}}\{2q{x}_0^2\left(D_{1\alpha }+{\displaystyle \frac{M_2}{x_0^3}}\right)+x_0\left(D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-zd)\right)\hfill \\ +gh\left(\frac{x_0^2}{3}+\frac{3x_0}{2\alpha }+\frac{2}{{\alpha }^2}+h\left(\frac{x_0}{12}+\frac{1}{6\alpha }\right)\right)\},\hfill \end{array}$$

it follows that

$$\begin{array}{c}{q}_{h{\alpha }_{op}}={\displaystyle \frac{-gh\left(\frac{x_0^2}{3}+\frac{3x_0}{2\alpha }+\frac{2}{{\alpha }^2}\right)-g{h}^2\left(\frac{x_0}{12}+\frac{1}{6\alpha }\right)-x_0\left(D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-zd)\right)}{2x_0^2\left(D_{1\alpha }+\frac{M_2}{x_0^3}\right)}}.\hfill \end{array}$$

Working algebraically we get formula (67).

b)

The estimate (68) follows straightforwardly from (67) with

$$C_{8\alpha }=\left|{\displaystyle \frac{g{B}_{1\alpha }}{6}}\right|.$$

From formula (65) we obtain that

$$\begin{array}{c}{J}_{2h\alpha }\left(q_{h{\alpha }_{op}}\right)-J_{2\alpha }\left(q_{{\alpha }_{op}}\right)=J_{2\alpha }\left(q_{h{\alpha }_{op}}\right)-J_{2\alpha }\left(q_{{\alpha }_{op}}\right)\hfill \\ +\frac{1}{2}{x}_0{y}_{0}gh(\frac{q_{h{\alpha }_{op}}{x}_0^2}{3}-\frac{5g{x}_0^3}{24}-\frac{(b-z_d)x_0}{2}+\frac{3x_0{q}_{h{\alpha }_{op}}}{2\alpha }\hfill \\ -\frac{7g{x}_0^2}{6\alpha }-\frac{2(b-z_d)}{\alpha }+\frac{2(q_{h{\alpha }_{op}}-g{x}_0)}{{\alpha }^2})+o\left(h^2\right).\hfill \end{array}$$

(70)

Moreover taking into account the explicit expression for $J_{2\alpha }$ given by (65) and formula (67) it follows that

$$\begin{array}{c}{J}_{2\alpha }\left(q_{h{\alpha }_{op}}\right)-J_{2\alpha }\left(q_{{\alpha }_{op}}\right)={\displaystyle \frac{x_0^2{y}_0}{2}}\left(q_{h{\alpha }_{op}}-q_{{\alpha }_{op}}\right)[\left(q_{h{\alpha }_{op}}+q_{{\alpha }_{op}}\right)x_0\left(D_{1\alpha }+\frac{M_2}{x_0^3}\right)\hfill \\ +D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-z_d)]\hfill \\ =-{\displaystyle \frac{x_0^2{y}_0}{12}}g{B}_{1\alpha }\left(2q_1{x}_0\left(D_{1\alpha }+{\displaystyle \frac{M_2}{x_0^3}}\right)+D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-z_d)\right)+o\left(h^2\right).\hfill \end{array}$$

(71)

Combining (70) and (71) we get estimate (69) with

$$\begin{array}{cc}{C}_{9\alpha }\hfill & ={\displaystyle \frac{x_0^2{y}_0}{2}}\left|g\right||-{\displaystyle \frac{B_{1\alpha }}{6}}\left(2q_{{\alpha }_{op}}{x}_0\left(D_{1\alpha }+{\displaystyle \frac{M_2}{x_0^3}}\right)+D_{2\alpha }g{x}_0^2+D_{3\alpha }(b-z_d)\right)\hfill \\ \\ \hfill & +{\displaystyle \frac{q_{{\alpha }_{op}}{x}_0}{3}}-{\displaystyle \frac{5g{x}_0^2}{24}}-{\displaystyle \frac{(b-z_d)}{2}}+{\displaystyle \frac{3q_{{\alpha }_{op}}}{2\alpha }}-{\displaystyle \frac{7g{x}_0}{6\alpha }}-{\displaystyle \frac{2(b-z_d)}{\alpha x_0}}+{\displaystyle \frac{2(q_{{\alpha }_{op}}-g{x}_0)}{{\alpha }^2{x}_0}}|.\hfill \end{array}$$

□

Lemma 14. 

Let us consider $u_{q_{{\alpha }_{op}}}$ the solution of (1) -(3) for $q=q_{{\alpha }_{op}}$ and $u_{h{q}_{h{\alpha }_{op}}}$ the discrete solution given in (26) for $h>0$ and $q=q_{h{\alpha }_{op}}$. Then we have that:

$$\begin{array}{c}\mathrm{a}){\parallel u_{h\alpha q_{h{\alpha }_{op}}}-u_{\alpha q_{{\alpha }_{op}}}\parallel }_H\approx C_{10\alpha }h,\hfill \\ \\ \mathrm{b}){∥\frac{\partial u_{h\alpha q_{h{\alpha }_{op}}}}{\partial x}-\frac{\partial u_{\alpha q_{{\alpha }_{op}}}}{\partial x}∥}_H\approx C_{11\alpha }h,\hfill \end{array}$$

(72)

Proof. 

Similarly to what was done in Lemma 12 it is obtained that

$$\begin{array}{c}{C}_{10\alpha }=\left|g\right|x_0\hfill \\ \sqrt{\frac{x_0{y}_0}{108}}\sqrt{B_{1\alpha }^2\left(\frac{3}{{\alpha }^2{x}_0^2}+\frac{3}{\alpha x_0}+1\right)+9\left(\frac{12}{{\alpha }^2{x}_0^2}+\frac{6}{\alpha x_0}+1\right)-3B_{1\alpha }\left(\frac{12}{{\alpha }^2{x}_0^2}+\frac{9}{\alpha x_0}+2\right)}\hfill \end{array}$$

$$C_{11\alpha }={\displaystyle \frac{\left|g\right|}{6}}\sqrt{x_0{y}_0(12-6B_{1\alpha }+B_{1\alpha }^2)}.$$

□

Remark 7. 

The constants verify that $C_{i\alpha }\to C_i$, when $\alpha \to \infty ,$ for each $i=7,\dots ,11.$

Remark 8. 

The double convergence when $(h,\alpha )\to (0,+\infty )$ of the optimal control of the problem $\left(P_{2h\alpha }\right)$ holds. The relationship of the optimal control problems $\left(P_2\right)$, $\left(P_{2\alpha }\right)$, $\left(P_{2h}\right)$ and $\left(P_{2h\alpha }\right)$ is given by the following diagram:

5. Boundary Optimization Problem with Variable b

5.1. Discrete Problem $\left(P_{3h}\right)$ Associated to $\left(P_3\right)$

In this section we consider the boundary optimal control problem $\left(P_3\right)$ given by (10). Taking into account expression (12), for a given constant b we get that

$$\begin{array}{c}{J}_3\left(b\right)=\frac{x_0{y}_0}{2}\{b^2\left(1+\frac{M_3}{x_0}\right)+b\left(\frac{2g{x}_0^2}{3}-q{x}_0-2z_d\right)\hfill \\ +\left[\frac{2g^2{x}_0^4}{15}-\frac{5gq{x}_0^3}{12}+\frac{x_0^2}{3}(q^2-2g{z}_d)+z_d(z_d+q{x}_0)\right]\}.\hfill \end{array}$$

(73)

Then the boundary optimal variable of the problem $\left(P_3\right)$, called $b_{op}$, and the associated continuous optimal state, are given respectively by:

$$b_{op}={\displaystyle \frac{-\frac{g{x}_0^2}{3}+\frac{q{x}_0}{2}+z_d}{1+\frac{M}{x_0}}},u_{b_{op}}=-\frac{1}{2}g{x}^2+(g{x}_0-q)x+b_{op}.$$

(74)

We define the discrete optimal control problem $\left(P_{3h}\right)$ on the constant temperature b as

$$\mathrm{find}{b}_{h_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{3h}\left(b_{h_{op}}\right)=\underset{b\in \mathbb{R}}{min}{J}_{3h}\left(b\right)$$

where the discrete cost function $J_{3h}\left(b\right)$ is defined as:

$$J_{3h}\left(b\right)=\frac{1}{2}{\parallel u_h-z_d\parallel }_H^2+\frac{1}{2}{M}_3{\parallel b\parallel }_B^2,$$

where $u_h$ is given in (22), h is the spatial step, and $z_d$ (the desired state) is constant.

Notice that the cost function $J_{3h}$ can be explicitly written as:

$$\begin{array}{c}{J}_{3h}\left(b\right)=J_3\left(b\right)+\frac{x_0{y}_{0}g}{2}\{-bh\left(\frac{x_0}{2}+\frac{h}{6}\right)+h{x}_0\left(\frac{1}{3}q{x}_0-\frac{5}{24}g{x}_0^2+\frac{z_d}{2}\right)\hfill \\ +\frac{1}{6}{h}^2\left(\frac{q{x}_0}{2}+\frac{1}{6}g{x}_0^2+z_d\right)+\frac{1}{24}g{h}^3{x}_0+\frac{1}{180}g{h}^4\}.\hfill \end{array}$$

(75)

Lemma 15. 

Let $b\in R$ and $h>0$ we have that:

$$|J_{3h}\left(b\right)-J_3\left(b\right)|\approx C_{12}h,$$

(76)

where

$$C_{12}=\frac{1}{2}{x}_0^2{y}_0\left|g\right|\left|-\frac{b}{2}+\frac{q{x}_0}{3}-\frac{5g{x}_0^2}{24}+\frac{z_d}{2}\right|,$$

does not depend on $h.$

Proof. 

It follows from expression (75) for $J_{3h}$. □

Lemma 16. 

Let us consider $h>0$.

a) 

The explicit expression for the optimal variable $b_{h_{op}}$ is given by:

$$b_{h_{op}}=b_{op}+E_{1}g{x}_{0}h\left(1+\frac{h}{3x_0}\right),E_1={\displaystyle \frac{1}{4\left(1+\frac{M_3}{x_0}\right)}}.$$

(77)

b) 

The following error estimates hold:

$$|b_{h_{op}}-b_{op}|\approx C_{13}h,$$

(78)

$$\left|J_{3h}\left(b_{h_{op}}\right)-J_3\left(b_{op}\right)\right|\approx C_{14}h,$$

(79)

where $C_{13}$ and $C_{14}$ do not depend on h.

Proof. 

a)

According to (75) we have that

$$J_{3h}^{\prime }\left(b\right)=J_3^{\prime }\left(b\right)-\frac{x_0{y}_0}{2}gbh\left(\frac{x_0}{2}+\frac{h}{6}\right).$$

(80)

Then, formula (77) for $b_{h_{op}}$ follows immediately.

b)

Estimate (78) is a direct consequence of (77) with $C_{13}=\left|E_{1}g{x}_0\right|$.

Moreover, taking into account formulas (73) and (75) for $J_3$ and $J_{3h}$ and formulas (74) and (77) for $b_{op}$ and $b_{h_{op}}$, respectively, it follows that

$$\begin{array}{c}{J}_{3h}\left(b_{h_{op}}\right)-J_3\left(b_{op}\right)\hfill \\ =J_3\left(b_{h_{op}}\right)-J_3\left(b_{op}\right)+\frac{x_0^2{y}_{0}g}{2}\left\{-\frac{b_{h_{op}}}{2}+\frac{1}{3}q{x}_0-\frac{5}{24}g{x}_0^2+\frac{z_d}{2}\right\}h+o\left(h^2\right).\hfill \end{array}$$

In addition, the expression $J_3\left(b_{h_{op}}\right)-J_3\left(b_{op}\right)$ can be rewritten as

$$J_3\left(b_{h_{op}}\right)-J_3\left(b_{op}\right)=\frac{x_0^2{y}_{0}g}{2}{E}_1\left(2b_{op}\left(1+\frac{M_3}{x_0}\right)+\frac{2}{3}g{x}_0^2-q{x}_0-2z_d\right)h+o\left(h^2\right).$$

Therefore it follows estimate (2) with

$$C_{14}=\frac{x_0^2{y}_{0}g}{2}\left|E_1\left(2b_{op}\left(1+\frac{M_3}{x_0}\right)+\frac{2}{3}g{x}_0^2-q{x}_0-2z_d\right)-\frac{b_{op}}{2}+\frac{q{x}_0}{3}-\frac{5g{x}_0^2}{24}+\frac{z_d}{2}\right|.$$

□

Lemma 17. 

Let us consider $u_{b_{op}}$ the solution of (1) -(3) for $b=b_{op}$ and $u_{h{b}_{h_{op}}}$ the discrete solution given in (26) for $h>0$ and $b=b_{h_{op}}.$ Then we have that:

$$\mathrm{a}){\parallel u_{h{b}_{h_{op}}}-u_{b_{op}}\parallel }_H\approx C_{15}h,\mathrm{b}){∥\frac{\partial u_{h{b}_{h_{op}}}}{\partial x}-\frac{\partial u_{b_{op}}}{\partial x}∥}_H\approx C_{16}h,$$

(81)

where $C_{15}$ and $C_{16}$ are constants that do not depend on h.

Proof. 

Working algebraically it is obtained that

$${\parallel u_{h{b}_{h_{op}}}-u_{b_{op}}\parallel }_H^2=\frac{x_0^3{y}_0{g}^2}{2}\left(2E_1^2-E_1+\frac{1}{6}\right)h^2+o\left(h^3\right).$$

Then it is obtained estimate $\mathrm{a})$ with

$$C_{15}=x_0\left|g\right|\sqrt{\frac{x_0{y}_0}{2}}\sqrt{2E_1^2-E_1+\frac{1}{6}}.$$

In a similar manner we get that estimate $\mathrm{b})$ holds with

$$C_{16}=\left|g\right|\sqrt{\frac{x_0{y}_0}{2}}.$$

□

5.2. Discrete Problem $\left(P_{3h\alpha }\right)$ Associated to $\left(P_{3\alpha }\right)$

From [6], we know that the continuous quadratic functional cost in (6) for the optimization problem $\left(P_{3\alpha }\right)$ is explicitly given by:

$$J_{3\alpha }\left(b\right)=J_3\left(b\right)+\frac{x_0{y}_0}{2}\left\{\frac{1}{3\alpha }(q-g{x}_0)(-6b+3q{x}_0-2g{x}_0^2+6z_d)+\frac{1}{{\alpha }^2}{(q-g{x}_0)}^2\right\}$$

(82)

where $J_3$ is defined by (73). Moreover, the continuous optimal boundary control $b_{{\alpha }_{op}}$ is given by

$$b_{{\alpha }_{op}}=b_{op}-{\displaystyle \frac{g{x}_0-q}{\alpha (1+\frac{M_3}{x_0})}}.$$

(83)

The continuous associated state is established by:

$$u_{b_{{\alpha }_{op}}}(x,y)=-\frac{1}{2}g{x}^2+(g{x}_0-q)x+\frac{1}{\alpha }(g{x}_0-q)+b_{{\alpha }_{op}}.$$

(84)

Defining the discrete cost function as:

$$J_{3h\alpha }\left(b\right)=\frac{1}{2}{\parallel u_{h\alpha }-z_d\parallel }_H^2+\frac{1}{2}{M}_3{\parallel b\parallel }_B^2$$

(85)

where $u_{h\alpha }$ is the solution of $\left(S_{h\alpha }\right)$ given in (26), we set the following discrete optimization problem $\left(P_{3h\alpha }\right)$ as

$$\mathrm{find}{b}_{h{\alpha }_{op}}\in \mathbb{R}\mathrm{such}\mathrm{that}{J}_{3h\alpha }\left(b_{h{\alpha }_{op}}\right)=\underset{b\in \mathbb{R}}{min}{J}_{3h\alpha }\left(b\right).$$

Working algebraically lead us to write $J_{3h\alpha }$ as follows:

$$\begin{array}{c}{J}_{3h\alpha }\left(b\right)=J_{3\alpha }\left(b\right)+\frac{1}{2}{x}_0{y}_{0}gh\{-b\left(\frac{x_0}{2}+\frac{2}{\alpha }+\frac{h}{6}\right)+g\left(\frac{-5x_0^3}{24}-\frac{7x_0^2}{6\alpha }-\frac{2x_0}{{\alpha }^2}\right)\hfill \\ +q\left(\frac{x_0^2}{3}+\frac{3x_0}{2\alpha }+\frac{2}{{\alpha }^2}\right)+z_d\left(\frac{x_0}{2}+\frac{2}{\alpha }\right)\hfill \\ +h\left[g\left(\frac{x_0^2}{36}+\frac{x_0}{3\alpha }+\frac{1}{{\alpha }^2}\right)+q\left(\frac{x_0}{12}+\frac{1}{6\alpha }\right)+\frac{z_d}{6}\right]+h^{2}g\left(\frac{x_0}{24}+\frac{1}{6\alpha }\right)+h^3\frac{g}{180}\}.\hfill \end{array}$$

(86)

Lemma 18. 

For $b\in B$ and $h>0$ we have that

$$|J_{3h\alpha }\left(b\right)-J_{3\alpha }\left(b\right)|\approx C_{12\alpha }h,$$

(87)

with

$$C_{12\alpha }=\frac{x_0{y}_0}{2}\left|g\right|\left|-b\left(\frac{x_0}{2}+\frac{2}{\alpha }+\frac{h}{6}\right)+g\left(\frac{-5x_0^3}{24}-\frac{7x_0^2}{6\alpha }-\frac{2x_0}{{\alpha }^2}\right)+q\left(\frac{x_0^2}{3}+\frac{3x_0}{2\alpha }+\frac{2}{{\alpha }^2}\right)+z_d\left(\frac{x_0}{2}+\frac{2}{\alpha }\right)\right|.$$

Proof. 

It arises immediately from (86). □

Lemma 19. 

Let us consider $h>0$.

a) 

The explicit expression for the optimal control $b_{h{\alpha }_{op}}$ is given by:

$$b_{h{\alpha }_{op}}=b_{{\alpha }_{op}}+E_{1}g{x}_{0}h\left(1+\frac{4}{\alpha x_0}+\frac{h}{3x_0}\right),$$

(88)

where $E_1$ is given in (77).

b) 

The following error estimates hold:

$$|b_{h{\alpha }_{op}}-b_{{\alpha }_{op}}|\approx C_{13\alpha }h,$$

(89)

$$\left|J_{3h\alpha }\left(b_{h{\alpha }_{op}}\right)-J_{3\alpha }\left(b_{{\alpha }_{op}}\right)\right|\approx C_{14\alpha }h,$$

(90)

where $C_{13\alpha }$ and $C_{14\alpha }$ do not depend on h.

Proof. 

a)

It follows from the expression $J_{3h\alpha }$ given by (86).

b)

The estimate in (89) is obtained immediately from item $\mathrm{a})$ with

$$C_{13\alpha }=|E_1\left|\right|g|x_0\left|1+\frac{4}{\alpha x_0}\right|.$$

Taking into account (86) and (82) it results that

$$J_{3h\alpha }\left(b_{h{\alpha }_{op}}\right)-J_{3\alpha }\left(b_{{\alpha }_{op}}\right)=J_{3\alpha }\left(b_{h{\alpha }_{op}}\right)-J_{3\alpha }\left(b_{{\alpha }_{op}}\right)+F_{1\alpha }h+o\left(h^2\right).$$

with

$$F_{1\alpha }=\frac{1}{2}{x}_0^2{y}_{0}g\left(-(b_{{\alpha }_{op}}-z_d)\left(\frac{1}{2}+\frac{2}{\alpha x_0}\right)+g{x}_0^2\left(\frac{-5}{24}-\frac{7}{6\alpha x_0}-\frac{2}{{\alpha }^2{x}_0^2}\right)+q{x}_0\left(\frac{1}{3}+\frac{3}{2\alpha x_0}+\frac{2}{{\alpha }^2{x}_0^2}\right)\right).$$

In addition, from the definition of $J_{3\alpha }$ and $b_{h{\alpha }_{op}}$ we have that

$$J_{3\alpha }\left(b_{h{\alpha }_{op}}\right)-J_{3\alpha }\left(b_{{\alpha }_{op}}\right)=J_3\left(b_{h{\alpha }_{op}}\right)-J_3\left(b_{{\alpha }_{op}}\right)+F_{2\alpha }h+o\left(h^2\right)$$

with

$$F_{2\alpha }=\frac{x_0^2{y}_{0}g}{\alpha }{E}_1(-q+g{x}_0)\left(1+\frac{4}{\alpha x_0}\right).$$

Finally, according to formula (73) for $J_3$ we get that

$$J_3\left(b_{h{\alpha }_{op}}\right)-J_3\left(b_{{\alpha }_{op}}\right)=F_{3\alpha }h+o\left(h^2\right),$$

with

$$F_{3\alpha }=\frac{x_0^2{y}_{0}g}{2}{E}_1\left(1+\frac{4}{\alpha x_0}\right)\left(2b_{{\alpha }_{op}}\left(1+\frac{M_3}{x_0}\right)+\frac{2}{3}g{x}_0^2-q{x}_0-2z_d\right).$$

Therefore, estimate (2) holds for

$$C_{14\alpha }=\left|F_{1\alpha }+F_{2\alpha }+F_{3\alpha }\right|.$$

□

Lemma 20. 

Let us consider $u_{b_{{\alpha }_{op}}}$ the solution of (1) -(3) for $b=b_{{\alpha }_{op}}$ and $u_{h{b}_{h{\alpha }_{op}}}$ the discrete solution given in (26) for $h>0$ and $b=b_{h{\alpha }_{op}}.$ Then we have that:

$$\mathrm{a}){\parallel u_{h\alpha b_{h{\alpha }_{op}}}-u_{\alpha b_{{\alpha }_{op}}}\parallel }_H\approx C_{15\alpha }h,\mathrm{b}){∥\frac{\partial u_{h\alpha b_{h{\alpha }_{op}}}}{\partial x}-\frac{\partial u_{\alpha b_{{\alpha }_{op}}}}{\partial x}∥}_H\approx C_{16\alpha }h,$$

(91)

Proof. 

Similarly to what was done in Lemma 12 it is obtained that

$$\begin{array}{c}{C}_{15\alpha }=x_0\left|g\right|\sqrt{\frac{x_0{y}_0}{2}}\mathcal{A},\\ \mathcal{A}=\sqrt{E_1^2\left(2+\frac{16}{\alpha x_0}+\frac{32}{{\alpha }^2{x}_0^2}\right)+E_1\left(-1-\frac{8}{\alpha x_0}-\frac{16}{{\alpha }^2{x}_0^2}\right)+\frac{1}{6}+\frac{1}{\alpha x_0}+\frac{2}{{\alpha }^2{x}_0^2}},\\ C_{16\alpha }=C_{16}.\end{array}$$

□

Remark 9. 

The constants obtained in the estimates of the previous lemmas verify that $C_{i\alpha }\to C_i$, when $\alpha \to \infty$ for $i=12,\dots ,16$.

Remark 10. 

The double convergence when $(h,\alpha )\to (0,+\infty )$ of the optimal control of the problem $\left(P_{3h\alpha }\right)$ holds. The relationship of the optimal control of the problems $\left(P_3\right)$, $\left(P_{3\alpha }\right)$, $\left(P_{3h}\right)$ and $\left(P_{3h\alpha }\right)$ is given by the following diagram:

6. Numerical Results

We will carry out some numerical simulations in order to illustrate the theoretical results obtained in the previous sections for the optimal control problems $\left(P_{ih}\right)$ and $\left(P_{ih\alpha }\right)$ for $i=1,2,3$.

Throughout this section we consider the domain $\Omega =[0,1]\times [0,1]$, i.e, $x_0=y_0=1$.

Before analyzing the optimal control problems we illustrate the behavior of the continuous state of the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$ and the discrete state of the systems $\left(S_h\right)$ and $\left(S_{h\alpha }\right)$.

In Figure 1 a) we plot the state of the system u given by (13) and the approximate discrete function $u_h$ defined by (22) against the position x for $h=1/3,1/5,1/10$. As we saw in Lemma 1 for each fixed x, the $u_h\left(x\right)$ increase and get closer to the limit $u\left(x\right)$ as h decreases. In a similar manner, in Figure 1 b), for $\alpha =50$ we obtain the system $u_{\alpha }$ given by (13) and the approximate discrete function $u_{h\alpha }$ defined by (26) against the position x for $h=1/3,1/5,1/10$. Notice that as h decreases, the functions $\{u_{h\alpha }\}$ increase and get closer to the limit $u_{\alpha }$ as it was proved in Lemma 2.

In addition in order to visualize the double convergence of $u_{h\alpha }\to u$ when $(h,\alpha )\to (0,\infty )$, in Figure 2 we plot u and $u_{h\alpha }$ for $(h,\alpha )=(1/3,10),(1/5,50)$ and $(1/10,500)$.

6.1. Control Variable g

In this subsection we obtain some computational examples for the optimal distributed control problems $\left(P_1\right)$, $\left(P_{1h}\right)$, $\left(P_{1\alpha }\right)$ and $\left(P_{1h\alpha }\right)$. For each plot we set $q=12,b=30,z_d=40$ and $M_1=1$.

In Figure 3 we plot the continuous quadratic cost function $J_1$ given by (27) and the discrete cost function $J_{1h}$ obtained in (32) against g for $h=1/10$, $1/50$ and $1/100$. Notice that as h decreases the function $J_{1h}=J_{1h}\left(g\right)$ also decreases to the limit function $J_1=J_1\left(g\right)$ in agreement with Lemma 3. In a similar manner in Figure 3, for $\alpha =50$ we obtain the continuous function $J_{1\alpha }$ and the discrete functions $J_{1h\alpha }$ for $h=1/10,1/50$ and $1/100$ observing the convergence of $J_{1h\alpha }\to J_{1\alpha }$ as h decreases to 0. Moreover, Figure 3 shows the double convergence of $J_{1h\alpha }\to J_1$ when $(h,\alpha )\to (0,\infty )$. We illustrate how $J_{1h\alpha }$ gets closer to $J_1$ as the value of h decreases and the value of $\alpha$ increases.

In Figure 3 we plot the continuous optimal control $g_{op}$ for the problem $\left(P_1\right)$ given by (29) and the optimal control $g_{{\alpha }_{op}}$ given by (41) for $\alpha =15,50,100$. Notice that as $\alpha$ increases, $g_{{\alpha }_{op}}$ decreases to the limit $g_{op}$. In addition, we set different values of n between $n=10$ and $n=100$. Recalling that $h=\frac{x_0}{n}=\frac{1}{n}$, for each h we obtain the optimal discrete control $g_{h_{op}}$ to the problem $\left(P_{1h}\right)$ defined by (4) and the optimal discrete control $g_{h{\alpha }_{op}}$ to the problem $\left(P_{1h\alpha }\right)$ given by (41) for $\alpha =15,50,100$. For each $\alpha$ fixed, we have that the discrete solution $g_{h{\alpha }_{op}}\to g_{{\alpha }_{op}}$ when $h\to 0$, i.e. $n\to \infty$.

6.2. Control Variable q

In this subsection we run some computational examples for the optimal boundary control problems $\left(P_2\right)$, $\left(P_{2h}\right)$, $\left(P_{2\alpha }\right)$ and $\left(P_{2h\alpha }\right)$. For each plot we set $g=10,b=50,z_d=40$ and $M_2=1$.

In Figure 4 we plot the continuous quadratic cost function $J_2$ given by (50) and the discrete cost function $J_{2h}$ obtained in (52) against q for $h=1/10$, $1/25$ and $1/50$. Observe that as h decreases the function $J_{2h}=J_{2h}\left(q\right)$ also decreases to the limit function $J_2=J_2\left(q\right)$. In a similar way, in Figure 4, for $\alpha =100$ we obtain the continuous function $J_{2\alpha }$ and the discrete functions $J_{2h\alpha }$ for $h=1/10,1/25$ and $1/50$. The convergences $J_{2h}\to J_2$ and $J_{2h\alpha }\to J_{2\alpha }$ when $h\to 0$ are in agreement with Lemmas 9 and 12, respectively.

Moreover, Figure 4 shows the double convergence of $J_{2h\alpha }\to J_2$ when $(h,\alpha )\to (0,\infty )$. We illustrate how $J_{2h\alpha }$ gets closer to $J_2$ as the value of h decreases and the value of $\alpha$ increases.

In Figure 4 we plot the continuous optimal control $q_{op}$ for the problem $\left(P_2\right)$ given by (51) and the optimal control $q_{{\alpha }_{op}}$ given by (63) for $\alpha =50,100,200$. Notice that as $\alpha$ increases, $q_{{\alpha }_{op}}$ decreases to the limit $q_{op}$. In addition, we set different values of n between $n=10$ and $n=100$. Recalling that $h=\frac{x_0}{n}=\frac{1}{n}$, for each h we obtain the optimal discrete control $q_{h_{op}}$ to the problem $\left(P_{2h}\right)$ defined by (54) and the optimal discrete control $q_{h{\alpha }_{op}}$ to the problem $\left(P_{2h\alpha }\right)$ given by (67) for $\alpha =50,100,200$. For each $\alpha$ fixed, we have that the discrete solution $q_{h{\alpha }_{op}}\to q_{{\alpha }_{op}}$ when $h\to 0$, i.e. $n\to \infty$.

6.3. Control Variable b

In this section we obtain some computational examples for the optimal distributed control problems $\left(P_3\right)$, $\left(P_{3h}\right)$, $\left(P_{3\alpha }\right)$ and $\left(P_{3h\alpha }\right)$. For each plot we set $q=12,g=10,z_d=40$ and $M_3=1$.

In Figure 5 we plot the continuous quadratic cost function $J_3$ given by (73) and the discrete cost function $J_{3h}$ obtained in (75) against g for $h=1/10$, $1/25$ and $1/100$. Notice that as h decreases the function $J_{3h}=J_{3h}\left(b\right)$ also decreases to the limit function $J_3=J_3\left(b\right)$ in agreement with Lemma 15. In a similar manner in Figure 5, for $\alpha =50$ we obtain the continuous function $J_{3\alpha }$ and the discrete functions $J_{3h\alpha }$ for $h=1/10,1/25$ and $1/100$. Observe the convergence of $J_{3h\alpha }\to J_{3\alpha }$ as $h\to 0$. Moreover, Figure 5 shows the double convergence of $J_{3h\alpha }\to J_3$ when $(h,\alpha )\to (0,\infty )$. We illustrate how $J_{3h\alpha }$ gets closer to $J_3$ as the value of h decreases and the value of $\alpha$ increases.

In Figure 5 we plot the continuous optimal control $b_{op}$ for the problem $\left(P_3\right)$ given by (74) and the optimal control $b_{{\alpha }_{op}}$ given by (83) for $\alpha =15,50,100$. Notice that as $\alpha$ increases, $b_{{\alpha }_{op}}$ decreases to the limit $b_{op}$. In addition, we set different values of n between $n=10$ and $n=100$. Recalling that $h=\frac{x_0}{n}=\frac{1}{n}$, for each h we obtain the optimal discrete control $b_{h_{op}}$ to the problem $\left(P_{3h}\right)$ defined by (77) and the optimal discrete control $b_{h{\alpha }_{op}}$ to the problem $\left(P_{3h\alpha }\right)$ given by (88) for $\alpha =15,50,100$. For each $\alpha$ fixed, we have that the discrete solution $b_{h{\alpha }_{op}}$ decreases to $b_{{\alpha }_{op}}$ when $h\to 0$.

7. Improvement of the Order of Convergence

In this section, we introduce alternative discrete solutions ${\tilde{u}}_h$ and ${\tilde{u}}_{h\alpha }$ associated with the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$, respectively, and analyze the order of convergence of ${\tilde{u}}_h$ to u and of ${\tilde{u}}_{h\alpha }$ to $u_{\alpha }$ as $h\to 0^{+}$. The Neumann boundary condition on ${\Gamma }_2$ is approximated by a three–point backward finite–difference scheme. Moreover, for the discrete solution ${\tilde{u}}_{h\alpha }$, the Robin boundary condition on ${\Gamma }_1$ is approximated by a three–point forward finite–difference scheme. These higher–order boundary approximations lead to an improved order of accuracy.

We consider the system $\left(S\right)$ defined by equations (1)–(2). From this system, we define the discrete problem $\left({\tilde{S}}_h\right)$, where ${\tilde{\mathrm{u}}}_i$ approximates $u(x_i,y)$, for $i=1,\dots ,n+1$. Notice that, from the Dirichlet condition on ${\Gamma }_1$ it follows immediately that ${\tilde{\mathrm{u}}}_1=b$.

For the interior nodes, we employ the classical centered second–order finite–difference approximation given in (15), which leads to the discrete system (16) for ${\tilde{\mathrm{u}}}_i$, $i=2,...,n$.

For the Neumann boundary condition on ${\Gamma }_2$, we use the three–point backward approximation

$$\frac{\partial u}{\partial x}(x_{n+1},y)\approx {\displaystyle \frac{3u(x_{n+1},y)-4u(x_n,y)+u(x_{n-1},y)}{2h}}.$$

(92)

Thus, the discrete Neumann condition can be written as

$$-2qh=3{\tilde{\mathrm{u}}}_{n+1}-4{\tilde{\mathrm{u}}}_n+{\tilde{\mathrm{u}}}_{n-1}.$$

(93)

In addition, from (16) for $i=n$, we obtain

$$-g{h}^2={\tilde{\mathrm{u}}}_{n+1}-2{\tilde{\mathrm{u}}}_n+{\tilde{\mathrm{u}}}_{n-1}.$$

(94)

Subtracting the two previous equations, it follows that

$$-{\tilde{\mathrm{u}}}_n+{\tilde{\mathrm{u}}}_{n+1}={\displaystyle \frac{g{h}^2}{2}}-qh.$$

(95)

Therefore the system given by (16) together with (95) can be written as

$$A\overline{\mathrm{w}}=T^{*}$$

(96)

where $\overline{\mathrm{w}}={\left({\tilde{\mathrm{u}}}_i\right)}_{i=2,\dots ,n+1}\in {\mathbb{R}}^n$ is the vector of unknowns, A is the matrix given by (20) and $T^{*}\in {\mathbb{R}}^n$ is the vector of independent terms:

$$T^{*}={\left(-g{h}^2-b,-g{h}^2,\dots ,-g{h}^2,-qh+\frac{g{h}^2}{2}\right)}^t.$$

(97)

Notice that the system (96) differs from (19) in the last component of the vector of independent terms. Solving the linear system gives

$${\tilde{\mathrm{u}}}_i=b+(i-1)h(g{x}_0-q)-{\displaystyle \frac{g{h}^2}{2}}{(i-1)}^2.$$

(98)

Taking into account that for $i=1,\dots ,n$

$${\tilde{m}}_i={\displaystyle \frac{{\tilde{\mathrm{u}}}_{i+1}-{\tilde{\mathrm{u}}}_i}{x_{i+1}-x_i}}=g{x}_0-q-(2i-1){\displaystyle \frac{gh}{2}},$$

(99)

and

$${\tilde{h}}_i={\tilde{\mathrm{u}}}_i-{\tilde{m}}_i{x}_i=b+{\displaystyle \frac{g{h}^2}{2}}i(i-1),$$

(100)

the linear approximation is given by ${\tilde{u}}_h(x,y)={\tilde{m}}_{i}x+{\tilde{h}}_i$, i.e.

$${\tilde{u}}_h(x,y)=\left(g{x}_0-q-(2i-1){\displaystyle \frac{gh}{2}}\right)x+{\displaystyle \frac{g{h}^2}{2}}i(i-1)+b,x\in [x_i,x_{i+1}],i=1,\dots ,n$$

(101)

In the following lemma, we give some bounds for the approximate function ${\tilde{u}}_h:$

Lemma 21. 

The following bounds hold:

$${\parallel u-{\tilde{u}}_h\parallel }_H\le D_1{h}^2,\mathrm{and}{\parallel \frac{\partial u}{\partial x}-\frac{\partial {\tilde{u}}_h}{\partial x}\parallel }_H\le {\tilde{D}}_{1}h,$$

where $D_1=\sqrt{{\displaystyle \frac{x_0{y}_0}{120}}}g$ and ${\tilde{D}}_1=\sqrt{{\displaystyle \frac{x_0{y}_0}{12}}}g$.

Proof. 

From the definition of the norm in the space H and using the expressions (13) and (101) for the functions u and ${\tilde{u}}_h$, respectively, it follows that

$$\begin{array}{cc}\hfill \parallel u-{\tilde{u}}_h{\parallel }_H^2& ={\int }_0^{y_0}{\int }_0^{x_0}{\left(u(x,y)-{\tilde{u}}_h(x,y)\right)}^{2}dxdy\hfill \\ \hfill & =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{E}_i^2\left(x\right)dx,\hfill \end{array}$$

(102)

where

$$E_i\left(x\right)=u(x,y)-{\tilde{u}}_h(x,y),x\in [x_i,x_{i+1}],y\in [0,y_0].$$

Note that, within each subinterval, $E_i\left(x\right)$ depends only on x and the index i, but not on y, since both u and ${\tilde{u}}_h$ are constant along the y-direction.

A direct computation yields

$$E_i\left(x\right)={\displaystyle \frac{g}{2}}\left(-x^2+(2i-1)hx-h^{2}i(i-1)\right)=-{\displaystyle \frac{g}{2}}(x-ih)(x-(i-1)h).$$

(103)

Then

$$\begin{array}{cc}\hfill {\int }_{x_i}^{x_{i+1}}{E}_i^2\left(x\right)dx& ={\displaystyle \frac{g^2}{4}}{\left[{\displaystyle \frac{{(x-ih)}^5}{5}}+{\displaystyle \frac{h}{2}}{(x-ih)}^4+{\displaystyle \frac{h^2}{3}}{(x-ih)}^3\right]}_{x_i}^{x_{i+1}}\hfill \\ \hfill & ={\displaystyle \frac{g^2}{4}}\left({\displaystyle \frac{h^5}{5}}-{\displaystyle \frac{h^5}{2}}+{\displaystyle \frac{h^5}{3}}\right)={\displaystyle \frac{g^2}{120}}{h}^5.\hfill \end{array}$$

(104)

As a consequence, from (102) it follows that

$$\parallel u-{\tilde{u}}_h{\parallel }_H^2=y_0{\displaystyle \frac{g^2{h}^5}{120}}n={\displaystyle \frac{x_0{y}_0{g}^2}{120}}{h}^4,$$

and then

$$\parallel u-{\tilde{u}}_h{\parallel }_H=\sqrt{{\displaystyle \frac{x_0{y}_0}{120}}}g{h}^2.$$

In addition,

$$\begin{array}{cc}\hfill \parallel \frac{\partial u}{\partial x}-\frac{\partial {\tilde{u}}_h}{\partial x}{\parallel }_H^2& ={\int }_0^{y_0}{\int }_0^{x_0}{\left(\frac{\partial u}{\partial x}(x,y)-\frac{\partial {\tilde{u}}_h}{\partial x}(x,y)\right)}^{2}dxdy\hfill \\ \hfill & =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{F}_i^2\left(x\right)dx,\hfill \end{array}$$

(105)

where

$$F_i\left(x\right)=\frac{\partial u}{\partial x}(x,y)-\frac{\partial {\tilde{u}}_h}{\partial x}(x,y)=-g\left(x-{\displaystyle \frac{(2i-1)h}{2}}\right),$$

for $x\in [x_i,x_{i+1}]$. Then

$$\begin{array}{cc}\hfill {\int }_{x_i}^{x_{i+1}}{F}_i^2\left(x\right)dx& ={\int }_{x_i}^{x_{i+1}}{g}^2{\left(x-{\displaystyle \frac{(2i-1)h}{2}}\right)}^{2}dx\hfill \\ & =g^2{\left[{\displaystyle \frac{1}{3}}{\left(x-{\displaystyle \frac{(2i-1)h}{2}}\right)}^3\right]}_{x_i}^{x_{i+1}}\hfill \\ & =g^2{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{h^3}{8}}+{\displaystyle \frac{h^3}{8}}\right)={\displaystyle \frac{g^2}{12}}{h}^3.\hfill \end{array}$$

Therefore, from (105) we have

$${∥\frac{\partial u}{\partial x}-\frac{\partial {\tilde{u}}_h}{\partial x}∥}_H^2=y_0{\displaystyle \frac{g^2}{12}}n{h}^3={\displaystyle \frac{x_0{y}_0{g}^2}{12}}{h}^2,$$

and finally

$${∥\frac{\partial u}{\partial x}-\frac{\partial {\tilde{u}}_h}{\partial x}∥}_H=\sqrt{{\displaystyle \frac{x_0{y}_0}{12}}}gh.$$

□

Remark 11. 

We emphasize that by improving the approximation of the Neumann boundary condition on ${\Gamma }_2$, the convergence order of the error $\parallel u-{\tilde{u}}_h{\parallel }_H$ is increased to second order, namely $O\left(h^2\right)$. This enhancement leads to a more accurate numerical approximation while remaining fully consistent with the theoretical convergence results established in [7,13].

Remark 12. 

The linear system (96) obtained by using the three–point backward finite–difference approximation for the Neumann boundary condition on ${\Gamma }_2$ can be equivalently interpreted by introducing a ghost point $x_{n+2}$ outside the computational domain and assuming that the discrete differential equation holds at the boundary node $x_{n+1}$. Indeed, assuming that the equation is satisfied at ${\tilde{\mathrm{u}}}_{n+1}$, we have

$$-g{h}^2={\tilde{\mathrm{u}}}_{n+2}-2{\tilde{\mathrm{u}}}_{n+1}+{\tilde{\mathrm{u}}}_n,$$

while the Neumann boundary condition is approximated by

$${\displaystyle \frac{{\tilde{\mathrm{u}}}_{n+2}-{\tilde{\mathrm{u}}}_n}{2h}}=-q.$$

Eliminating the ghost value ${\tilde{\mathrm{u}}}_{n+2}$ from these two expressions yields

$$-{\tilde{\mathrm{u}}}_n+{\tilde{\mathrm{u}}}_{n+1}=-qh+{\displaystyle \frac{g{h}^2}{2}},$$

which coincides with the boundary equation obtained in (95). Hence, the three–point backward finite–difference approximation of the Neumann condition is consistent with the ghost–point formulation and leads to the same discrete system.

Analogously to the analysis of system $\left(S\right)$, we propose a new discrete approximation ${\tilde{u}}_{h\alpha }$ for system $\left(S_{\alpha }\right)$ and study the order of convergence of ${\tilde{u}}_{h\alpha }$ to $u_{\alpha }$ as $h\to 0^{+}$. The associated discrete system $\left({\tilde{S}}_{h\alpha }\right)$ employs a three–point backward finite–difference approximation for the Neumann boundary condition on ${\Gamma }_2$ and a three–point forward finite–difference approximation for the Robin boundary condition on ${\Gamma }_1$, leading to improved accuracy.

We consider the system $\left(S_{\alpha }\right)$ defined by equations (1)–(3) and define ${\tilde{\mathrm{u}}}_{\alpha ,i}\approx u_{\alpha }(x_i,y)$.

For the interior nodes, $i=2,\dots ,n$, we employ the classical centered second–order finite–difference approximation given in (15):

$${\tilde{\mathrm{u}}}_{\alpha ,i+1}-2{\tilde{\mathrm{u}}}_{\alpha ,i}+{\tilde{\mathrm{u}}}_{\alpha ,i-1}=-g{h}^2.$$

(106)

For the Robin boundary at ${\Gamma }_1$, we use the three–point forward approximation:

$${\displaystyle \frac{-3{\tilde{\mathrm{u}}}_{\alpha ,1}+4{\tilde{\mathrm{u}}}_{\alpha ,2}-{\tilde{\mathrm{u}}}_{\alpha ,3}}{2h}}=\alpha ({\tilde{\mathrm{u}}}_{\alpha ,1}-b).$$

(107)

Combining this expression with the interior equation at $i=2$ yields the simplified discrete condition

$$-(1+\alpha h){\tilde{\mathrm{u}}}_{\alpha ,1}+{\tilde{\mathrm{u}}}_{\alpha ,2}=-\alpha hb-{\displaystyle \frac{g{h}^2}{2}}.$$

(108)

For the Neumann boundary at ${\Gamma }_2$ we use the three–point backward approximation:

$$3{\tilde{\mathrm{u}}}_{\alpha ,n+1}-4{\tilde{\mathrm{u}}}_{\alpha ,n}+{\tilde{\mathrm{u}}}_{\alpha ,n-1}=-2qh.$$

(109)

Combining with the interior equation for $i=n$ gives

$$-{\tilde{\mathrm{u}}}_{\alpha ,n}+{\tilde{\mathrm{u}}}_{\alpha ,n+1}=-qh+{\displaystyle \frac{g{h}^2}{2}}.$$

(110)

The system given by (106), (108) and (110) can be rewritten as

$$A_{\alpha }\overline{{\mathrm{w}}_{\alpha }}=T_{\alpha }^{*}$$

(111)

where ${\overline{\mathrm{w}}}_{\alpha }={\left({\widetilde{{\mathrm{u}}_{\alpha }}}_i\right)}_{i=1,\dots ,n+1}\in {\mathbb{R}}^{n+1}$ is the vector of unknowns, $A_{\alpha }$ is the matrix given by (24) and $T_{\alpha }^{*}\in {\mathbb{R}}^{n+1}$ is the vector of independent terms:

$$T_{\alpha }^{*}={\left(-\alpha bh-\frac{g{h}^2}{2},-g{h}^2,\dots ,-g{h}^2,-qh+\frac{g{h}^2}{2}\right)}^t\in {\mathbb{R}}^{n+1}.$$

(112)

It should be noted that only the first and last components of $T_{\alpha }^{*}$ differ from those in $T_{\alpha }$ given by (25).

The solution of the system (111) is given by

$${\tilde{u}}_{\alpha ,i}=b+{\displaystyle \frac{1}{\alpha }}\left(g{x}_0-q\right)+(i-1)h(g{x}_0-q)-{\displaystyle \frac{g}{2}}{\left((i-1)h\right)}^2,i=1,\dots ,n+1.$$

(113)

We define the linear interpolation on each subinterval $[x_i,x_{i+1}]$ by

$${\tilde{u}}_{h\alpha }(x,y)={\tilde{m}}_{\alpha ,i}x+{\tilde{h}}_{\alpha ,i},x\in [x_i,x_{i+1}],y\in [0,y_0],$$

(114)

where

$$\begin{array}{cc}\hfill {\tilde{m}}_{\alpha ,i}& =g{x}_0-q-gh\left(i-{\displaystyle \frac{1}{2}}\right),i=1,\dots ,n,\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill {\tilde{h}}_{\alpha ,i}& =b+{\displaystyle \frac{1}{\alpha }}(g{x}_0-q)+{\displaystyle \frac{g{h}^2}{2}}(i-1)i,i=1,\dots ,n.\hfill \end{array}$$

(116)

From the previous expressions, we derive the following lemma.

Lemma 22. 

The following bounds hold:

$${\parallel u_{\alpha }-{\tilde{u}}_{h\alpha }\parallel }_H\le D_2{h}^2,{\parallel \frac{\partial u_{\alpha }}{\partial x}-\frac{\partial {\tilde{u}}_{h\alpha }}{\partial x}\parallel }_H\le {\tilde{D}}_{2}h,$$

where $D_2=\sqrt{{\displaystyle \frac{x_0{y}_0}{120}}}g$ and ${\tilde{D}}_2=\sqrt{{\displaystyle \frac{x_0{y}_0}{12}}}g$.

Proof. 

By the definition of the H-norm, and using the expression for $u_{\alpha }$ in (13) as well as the definition of ${\tilde{u}}_{h\alpha }$ in (114), it follows that

$${\parallel u_{\alpha }-{\tilde{u}}_{h\alpha }\parallel }_H^2=y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{E}_{\alpha ,i}^2\left(x\right)dx,$$

(117)

where

$$E_{\alpha ,i}\left(x\right)=u_{\alpha }(x,y)-{\tilde{u}}_{h\alpha }(x,y)$$

$$E_{\alpha ,i}\left(x\right)=-{\displaystyle \frac{g}{2}}{x}^2+{\displaystyle \frac{g}{2}}\left(2i-1\right)hx-{\displaystyle \frac{g{h}^2}{2}}(i^2-i),x\in [x_i,x_{i+1}]$$

(118)

We can notice that $E_{\alpha ,i}\left(x\right)=E_i\left(x\right)$ where $E_i\left(x\right)$ is given by (103). Therefore, from (104), it follows immediately that

$${\parallel u_{\alpha }-{\tilde{u}}_{h\alpha }\parallel }_H^2\le {\displaystyle \frac{x_0{y}_0{g}^2}{120}}{h}^4.$$

In addition,

$$\begin{array}{cc}\hfill {∥\frac{\partial u}{\partial x}-\frac{\partial {\tilde{u}}_h}{\partial x}∥}_H^2& ={\int }_0^{y_0}{\int }_0^{x_0}{\left(\frac{\partial u}{\partial x}(x,y)-\frac{\partial {\tilde{u}}_h}{\partial x}(x,y)\right)}^{2}dxdy\hfill \\ \hfill & =y_0\sum _{i=1}^n{\int }_{x_i}^{x_{i+1}}{g}^2{\left(x-{\displaystyle \frac{(2i-1)h}{2}}\right)}^{2}dx\hfill \\ \hfill & =y_0{\displaystyle \frac{g^2}{12}}{h}^{3}n={\displaystyle \frac{x_0{y}_0{g}^2}{12}}{h}^2.\hfill \end{array}$$

(119)

□

8. Conclusions

Applying the finite difference method, we have derived the discrete systems $\left(S_h\right)$ and $\left(S_{h\alpha }\right)$ and the discrete optimization problems $\left(P_{ih}\right)$ and $\left(P_{ih\alpha }\right)$, $i=1,2,3$ where $\alpha >0$ is a parameter that represents the heat transfer coefficient on a portion of the boundary of the domain. Explicit discrete solutions have been found and convergence results when the discrete step h goes to zero and when $\alpha$ goes to infinite have been proved. Error estimations have been also obtained as a function of the step h. Some numerical computations have been provided in order to illustrate the theoretical results.

Finally, for the systems $\left(S\right)$ and $\left(S_{\alpha }\right)$, an alternative discretization of the Neumann boundary condition on ${\Gamma }_2$ and of the Robin boundary condition on ${\Gamma }_1$ for $\left(S_{\alpha }\right)$ has been considered. By modifying the approximation of these boundary conditions, the order of convergence of the numerical solution is improved, leading to a more accurate approximation.