Abstract

In this paper, we investigate a meshless approximation, the interpolating element-free Galerkin method, for an optimal control problem governed by fourth-order parabolic partial differential equations. The state, co-state and control variables are spatially discretized by an improved moving least squares approximation that satisfies the interpolation property, and time is discretized by a backward-Euler method. We derive some a priori error estimates for both the control and state approximations. Numerical experiments are presented to verify the theoretical results.

Keywords

Interpolating Element-Free Galerkin Method, Optimal Control Problem, Moving Least Squares Approximation, A Priori Error Estimates

Share and Cite:

1. Introduction

The optimal control problem is a topic problem in the field of computational mathematics, which has a wide range of applications in the fields of ecosystem, natural science, engineering technology, economic decision making, and so on. Based on Lions’s [1] systematic study of the optimal control problems for partial differential equations (PDEs), many scholars carried out research work on the mathematical theories and numerical methods, achieving various progress [2]-[4]. Meshless method, as an advanced numerical computation method, breaks the limitation of traditional mesh dependency [5] and solves PDEs based on approximation at discrete points. The meshless method is primarily categorized into two approaches: collocation-based and Galerkin-based. The Galerkin-based method has been the subject of extensive application and study due to its superior stability [6]. The element-free Galerkin (EFG) method, proposed by Belytschko et al. [7], evolves into the mainstream meshless methods, attributed to multiple factors, including its weak mesh dependency and exceptional stability. This method employs a global variational principle and utilizes a moving least squares (MLS) approximation to generate trial functions [7] [8].

There are numerous subsequent research developments based on the EFG method. Zhang et al. [9] proposed the improved moving least squares-Ritz meshless method, which substantially improved the computational efficiency while ensuring the computational accuracy. Dehghan’s team [10] successfully applied the interpolating element-free Galerkin (IEFG) method to solving the nonlinear Benjamin-Bona-Mahony-Burgers equation and the regular long-wave equation, thereby expanding the application scope of the method. Li and other academics [8] [11] applied the MLS approximation to propose a stable and improved algorithm, providing a theoretical guarantee of the reliability. Abbaszadeh et al. [12] innovatively combined the alternating direction implicit method with the IEFG method, successfully applying it to the solution of fractional order differential equations.

In recent years, meshless methods have made significant progress at both the theory and application. Fu’s research team [13] achieved tremendous success on the EFG method, and innovatively proposed a power analysis method based on the theory of regenerative smooth gradients. Meng et al. [14] successfully established a discrete scheme that can directly impose intrinsic boundary conditions by constructing a new type of approximation function with accurate interpolation property. Zhu and other academics [15] innovatively introduced the concept of spatio-temporal discretization into the radial basis meshless method, and proposed the spatio-temporal MQ radial basis method and a new algorithm coupled with polynomial basis for telegraph equations. Authors of [16] developed a specialized EFG method for obstacle problems in engineering, expanding the application scope of the method. Li and other scholars [17] conducted a systematic investigation into the algorithm for coupling the EFG method with the Bathe time-discrete scheme, which achieved monotonically convergent numerical solutions in elastic dynamics problems by optimizing the node arrangement and shape function support domains. Li et al. [18] [19] achieved significant results in theoretical analysis, rigorously proved the existence and uniqueness of the solution of the EFG method when the regenerative kernel gradient was used for smooth integrals, and established a complete theoretical system for the error estimates. In addressing the intricacies of quantum mechanics, Cui and other academics [20] developed an efficient EFG method, which was based on interpolating moving least squares (IMLS). The improved IEFG method proposed by the team of Wang [21], furnished a novel solution to the numerical simulation of elastic-plasticity problems by means of the introduction of a non-singular weight function to construct the shape function.

However, the IEFG method has not been widely adopted in the study of the optimal control problems governed by fourth-order parabolic PDEs. Duan [22] studied the optimal control problem of the extended Fisher-Kolmogorov (EFK) equation, which is a fourth-order nonlinear parabolic PDE. Although the scheme for optimal control problems under boundary conditions is given and the existence of the optimal solution to the equation is proven, no related research on numerical methods has been conducted.

Model Problem

In this paper, we study the following initial boundary problems for fourth-order parabolic PDEs. Let $\Omega \subset {ℝ}^d\left(d\ge 2\right)$ be a bounded region with a Lipschitz continous boundary $\partial \Omega$ , and let $J=\left[0,T\right]$ be the time interval. We seek a real-valued function $y=y\left(x,t\right)$ defined on $\Omega \times J$ , which satisfies

$y_t+{\kappa }_1{\Delta }^{2}y-{\kappa }_2\Delta y=f,x\in \partial \Omega ,t\in J,$ (1.1)

where ${\kappa }_1$ , ${\kappa }_2>0$ are given positive material parameters, $f$ is a given source term, and $\Delta$ denotes the Laplace operator. The operator ${\Delta }^2$ is the biharmonic operator, defined as the application of Laplace operator twice, i.e., ${\Delta }^{2}y=\Delta \left(\Delta y\right)$ . For the sake of simplicity, we define $f$ as the linear term.

Based on physical considerations, Equation (1.1) is supplemented with the following boundary condition:

$y\left(x,t\right)=\Delta y\left(x,t\right)=0,x\in \partial \Omega ,t\in J,$ (1.2)

and initial condition:

$y\left(x,0\right)=y_0\left(x\right),x\in \Omega .$ (1.3)

This paper establishes the optimal control problem based on fourth-order parabolic PDEs, and solves the equations by IEFG method which based on the IMLS approximation. The rest of this paper is organized as follows: in Section 2, the notations used throughout the paper and several theoretical results concerning the existence and uniqueness of solutions are introduced, and the optimality conditions are derived. In Section 3, we introduce the specific process to shape function construction of the IEFG method in detail, and establish the fully discrete approximation scheme for the continuous optimal control problem. In Section 4, the a priori error estimates for the control, state and co-state variables are provided under some specific assumptions. The results of the numerical experiments are presented in Section 5 in order to support the theoretical results. Finally, we summarize the conclusions.

2. Preliminaries and Optimality Conditions

In this section, we introduce some common notations, basic inequalities, and fundamental lemmas which will be utilized in the subsequent analysis. Then we proceed to a discussion of the optimal control problem governed by Equation (1.1), followed by an analysis of the existence and uniqueness of weak solutions to these equations. Ultimately, the optimality conditions are derived.

2.1. Notations

This subsection is intended to provide definitions for a series of symbols. For the sake of brevity, the symbols with equivalent meanings that are mentioned subsequently are not repeated.

Assume that $C$ is a positive constant which is independent of discrete parameters and may indicate different values in different circumstances. ${ℝ}^d$ is a $d$ -dimensional Euclidean space, and in this paper we consider the two-dimensional case (i.e., $d=2$ ).

Let $\Omega$ , ${\Omega }_U\subset {ℝ}^2$ be a bounded convex polygonal region with Lipschitz boundary $\partial \Omega$ and $\partial {\Omega }_U$ respectively, where $\Omega$ is the space to the state variable and ${\Omega }_U$ is the space to the control variable. Let $m$ be a non-negative integer, $1\le p\le \infty$ , and $W^{m,p}\left(\Omega \right)$ be a Sobolev space whose norm is denoted as ${\Vert \cdot \Vert }_{W^{m,p}\left(\Omega \right)}$ . To simplify the expression, $H^m\left(\Omega \right)$ is used to denote $W^{m,2}\left(\Omega \right)$ , with a norm shortened to ${\Vert \cdot \Vert }_m={\Vert \cdot \Vert }_{W^{m,2}\left(\Omega \right)}$ , and a semi-norm shortened to ${\left|\cdot \right|}_m={\left|\cdot \right|}_{W^{m,2}\left(\Omega \right)}$ . Define notation $H_0^m\left(\Omega \right)=\left\{v\in H^m\left(\Omega \right):{\frac{{\partial }^{s}v}{\partial x^s}|}_{\partial \Omega }=0,0\le s\le m\right\}$ . By convention, $\left(\cdot ,\cdot \right)$ is used to denote the $L^2$ -inner product, and it is briefly noted that $\Vert \cdot \Vert ={\Vert \cdot \Vert }_0$ is used to represent the $L^2$ -norm.

Denote $L^s\left(J;W^{m,p}\left(\Omega \right)\right)$ as the Banach space of all $L^s$ integrable functions from $J=\left[0,T\right]$ into $W^{m,p}\left(\Omega \right)$ with norm defined by

${\Vert \cdot \Vert }_{L^s\left(J;W^{m,p}\left(\Omega \right)\right)}={\left({\displaystyle {\int }_0^T{\Vert \cdot \Vert }_{W^{m,p}\left(\Omega \right)}^s\text{d}t}\right)}^{1/s},s\in \left[1,\infty \right),$

and with the standard modification

${\Vert \cdot \Vert }_{L^{\infty }\left(J;W^{m,p}\left(\Omega \right)\right)}=\text{ess}\underset{\left[0,T\right]}{\text{sup}}{\Vert \cdot \Vert }_{W^{m,p}\left(\Omega \right)},s=\infty .$

Let $0=t_0<t_1<\cdots <t_N=T$ be the division of time $J$ , and denote $k_i=t_i-t_{i-1} \left(i=1,\cdots ,N\right)$ , $k=\max k_i$ . The following notations are defined to make writing easier

$d_t{v}^i=\frac{v^i-v^{i-1}}{k_i},{\overline{d}}_t{v}^{i-1}=\frac{v^{i-1}-v^i}{k_i}.$(2.1)

We define the discrete time-dependent norms and the standard modification as follows:

${\left|\Vert v\Vert \right|}_{L^s\left(J;H^m\left(\Omega \right)\right)}:={\left({\displaystyle \sum }_{i=1}^N{k}_i{\Vert v^i\Vert }_m^s\right)}^{\frac{1}{s}},s=\left[1,\infty \right),$ (2.2)

${\left|\Vert v\Vert \right|}_{L^{\infty }\left(J;H^m\left(\Omega \right)\right)}:=\underset{i=1,\cdots ,N}{\text{sup}}{\Vert v^i\Vert }_m.$ (2.3)

2.2. Fundamental Lemma

Lemma 1. [23] Let ${\theta }_n$ , ${\varphi }_n$ and ${\delta }_n$ be non-negative sequences and $c_0>0$ , and assume that the sequence ${\theta }_n$ satisfies

$\{\begin{array}{l}{\theta }_0\le c_0,\\ {\theta }_n\le c_0+{\displaystyle \sum }_{k=0}^{n-1}{\delta }_k+{\displaystyle \sum }_{k=0}^{n-1}{\varphi }_k{\theta }_k,n\ge 1.\end{array}$

If $c_0\ge 0$ and ${\delta }_0\ge 0$ , then there are

${\theta }_n\le \left(c_0+{\displaystyle \sum }_{k=0}^{n-1}{\delta }_k\right)\exp \left({\displaystyle \sum }_{k=0}^{n-1}{\varphi }_k\right),n\ge 1.$

Here we briefly introduce two basic inequalities, which will be used repeatedly in the proof below.

$a\left(a-b\right)\ge \frac{1}{2}\left(a^2-b^2\right),\forall a,b\in ℝ,$ (2.4)

$\left|ab\right|\le \frac{\epsilon }{2}{a}^2+\frac{1}{2\epsilon }{b}^2,\forall a,b,\epsilon \in ℝ,\epsilon >0.$ (2.5)

2.3. Optimal Control Problems and Optimality Conditions

We now describe precisely the mathematical model of the control problems governed by Equation (1.1). To clarify the idea, we take the state space $W=L^2\left(0,T;V\right)$ , where $V=H^2\left(\Omega \right)\cap H_0^1\left(\Omega \right)$ , the control space $X=L^2\left(0,T;U\right)$ , where $U=L^2\left({\Omega }_U\right)$ , and the observation space $Y=L^2\left(0,T;L^2\left(\Omega \right)\right)$ . Let $B$ be a linear continuous operator from $X$ to $L^2\left(0,T;V^{\prime }\right)$ , and $B$ is the identity operator. Let $K$ be a closed convex set in $X$ :

$K=\left\{v\in X:v\ge 0,a.e.{\Omega }_U\times \left(0,T\right]\right\}.$

We focus on the following distributed optimal control problem (OCP):

$\underset{u\in K}{\text{min}}\left(\frac{1}{2}{\displaystyle {\int }_0^T\left({\Vert y-y_d\Vert }_{0,\Omega }^2+\alpha {\Vert u\Vert }_{0,{\Omega }_U}^2\right)\text{d}t}\right),$ (2.6)

$\{\begin{array}{ll}\frac{\partial y}{\partial t}+{\kappa }_1{\Delta }^{2}y-{\kappa }_2\Delta y=f+Bu,\hfill & x\in \Omega ,t\in \left(0,T\right],\hfill \\ y\left(x,t\right)=\Delta y\left(x,t\right)=0,\hfill & x\in \partial \Omega ,t\in \left[0,T\right],\hfill \\ y\left(x,0\right)=y_0\left(x\right),\hfill & x\in \Omega ,\hfill \end{array}$ (2.7)

where $\alpha$ is a positive number, $f\in L^2\left(0,T;L^2\left(\Omega \right)\right)$ , $y_d\in H^1\left(0,T;L^2\left(\Omega \right)\right)$ , and $y_0\in V$ .

Definition 2. [22] For all $w\in V$ , $t\in \left(0,T\right)$ , if the function $y\left(x,t\right)\in W$ is a weak solution to the problem (2.7). The problem is solved if it satisfies

$\left(\frac{\partial y}{\partial t},w\right)+{\kappa }_1\left(\Delta y,\Delta w\right)+{\kappa }_2\left(\nabla y,\nabla w\right)=\left(f+Bu,w\right).$

We present Theorem 3, which establishes the existence and uniqueness of weak solutions to (2.7).

Theorem 3. [22] Assuming that $y_0\in V$ then the problem (2.7) exists a unique weak solution $y\left(x,t\right)\in W$ .

Definitions of the inner products, bilinear forms, and related norms are given as follows:

$a\left(v,w\right)={\displaystyle {\int }_{\Omega }\nabla v\cdot \nabla w},{\Vert v\Vert }_a={\left({\displaystyle {\int }_{\Omega }\nabla v\cdot \nabla v}\right)}^{\frac{1}{2}},\forall v,w\in H^1\left(\Omega \right),$

$G\left(v,w\right)={\displaystyle {\int }_{\Omega }\Delta v\cdot \Delta w},{\Vert v\Vert }_g={\left({\displaystyle {\int }_{\Omega }\Delta v\cdot \Delta v}\right)}^{\frac{1}{2}},\forall v,w\in H^2\left(\Omega \right),$

$\left(v,w\right)={\displaystyle {\int }_{\Omega }v\cdot w},{\Vert v\Vert }_{0,\Omega }={\left({\displaystyle {\int }_{\Omega }{v}^2}\right)}^{\frac{1}{2}},\forall v,w\in L^2\left(\Omega \right),$

${\left(v,w\right)}_U={\displaystyle {\int }_{{\Omega }_U}v\cdot w},{\Vert v\Vert }_{0,{\Omega }_U}={\left({\displaystyle {\int }_{{\Omega }_U}{v}^2}\right)}^{\frac{1}{2}},\forall v,w\in L^2\left({\Omega }_U\right),$

${\left(v,w\right)}_V=a\left(v,w\right)+G\left(v,w\right),{\Vert v\Vert }_V={\left({\Vert v\Vert }_a^2+{\Vert v\Vert }_g^2\right)}^{\frac{1}{2}},\forall v,w\in V.$

To simplify the analysis, this paper takes ${\kappa }_1=1$ and ${\kappa }_2=1$ to determine Equation (1.1). The weak form of the state equation is as follows. For a given $f,u,y_0$ , find $y\left(u\right)\in H^1\left(0,T;L^2\left(\Omega \right)\right)\cap W$ for the (OCP) problem:

$\{\begin{array}{l}\left(\frac{\partial }{\partial t}y\left(u\right),w\right)+\left(\Delta y\left(u\right),\Delta w\right)+a\left(y\left(u\right),w\right)=\left(f+Bu,w\right),\\ \forall w\in V,t\in \left(0,T\right],\\ y\left(u\right)\left(x,0\right)=y_0\left(x\right),x\in \Omega ,\end{array}$ (2.8)

and we can solve the problem (2.8) uniquely.

Now introduce an objective functional:

$J\left(v\right)=\frac{1}{2}{\displaystyle {\int }_0^T\left({\Vert y\left(v\right)-y_d\Vert }_{0,\Omega }^2+\alpha {\Vert v\Vert }_{0,{\Omega }_U}^2\right)\text{d}t}.$

The above convex optimal control problem can be reformulated as follows, which labeled (QCP): find $u\in K$ such that

$J\left(u\right)=\underset{v\in K}{\min }J\left(v\right),$

which $y\left(v\right)\in W$ satisfies

$\{\begin{array}{l}\left(\frac{\partial }{\partial t}y\left(v\right),w\right)+\left(\Delta y\left(v\right),\Delta w\right)+a\left(y\left(v\right),w\right)=\left(f+Bv,w\right),\\ \forall w\in V,t\in \left(0,T\right],\\ y\left(v\right)\left(x,0\right)=y_0\left(x\right),x\in \Omega .\end{array}$

Based on Lions’s [1] theoretical analysis of optimal control problems, we obtain that the optimality condition for $u$ is given by the variational inequality

$J^{\prime }\left(u\right)\left(v-u\right)\ge 0,\forall v\in K,$

where $J^{\prime }\left(u\right)$ denotes the Gateaux derivative of $J\left(v\right)$ at $v=u$ . The following Lemma 4 is crucial in deriving the necssary optimality condition.

Lemma 4 [22] Suppose $f$ is a given state variable independent function, the mapping $v\to y\left(v\right)$ from $L^2\left(0,T;L^2\left({\Omega }_U\right)\right)$ to $W$ is weakly Gateaux-neutralizable at $v=u$ , such that the Gateaux derivative of $y\left(v\right)$ at $v=u$ in the direction $v-u\in L^2\left(0,T;L^2\left({\Omega }_U\right)\right)$ , denoted as $z=\mathcal{D}y\left(u\right)\left(v-u\right)$ , is the only weak solution to the following problem

$\{\begin{array}{ll}{z}_t+{\Delta }^{2}z-\Delta z=B\left(v-u\right),\hfill & \left(x,t\right)\in \Omega \times \left(0,T\right),\hfill \\ z\left(x,t\right)=\Delta z\left(x,t\right)=0,\hfill & \left(x,t\right)\in \partial \Omega \times \left(0,T\right),\hfill \\ z\left(0\right)=0,\hfill & x\in \Omega .\hfill \end{array}$

We utilize Lion’s optimal control theory to derive the following results. The control problem (QCP) has a unique solution $\left(y,u\right)$ , and $\left(y,u\right)$ is the solution of (QCP) if and only if there exists a co-state $p\in W$ such that the ternary $\left(y,p,u\right)$ satisfies the following optimality condition (QCP - OPT):

$\{\begin{array}{l}\left(\frac{\partial y}{\partial t},w\right)+\left(\Delta y,\Delta w\right)+a\left(y,w\right)=\left(f+Bu,w\right),\forall w\in V,\\ y\left(0\right)=y_0,\end{array}$ (2.9)

$\{\begin{array}{l}-\left(\frac{\partial p}{\partial t},q\right)+\left(\Delta p,\Delta q\right)+a\left(q,p\right)=\left(y-y_d,q\right),\forall q\in V,\\ p\left(T\right)=0,\end{array}$ (2.10)

${\displaystyle {\int }_0^T{\left(\alpha u+B^{*}p,v-u\right)}_U\text{d}t}\ge 0,\forall v\in K\subset X=L^2\left(0,T;U\right),$ (2.11)

where $B^{\text{*}}$ is the adjoint operator of $B$ and ${\left(\cdot ,\cdot \right)}_U$ is the inner product on ${\Omega }_U$ .

The inequality (2.11) is equivalent to

$u=\text{max}\left\{0,-\frac{1}{\alpha }{B}^{\text{*}}p\right\}.$ (2.12)

3. Meshless Basis Functions

3.1. The Interpolating MLS (IMLS) with Weight Function

Let $X=\left\{x_1,x_2,\cdots ,x_{N^{\text{*}}}\right\}$ be the set of all nodes in the bounded region $\Omega \subset {ℝ}^d$ where $N^{\text{*}}$ is the number of nodes. The parameter ${\rho }_I$ denotes the radius of the support domain of the node $x_I$ , and $\Vert \cdot \Vert$ denotes the Euclidean norm. The support domain of $x_I$ is defined as ${\Omega }_I=\left\{x|\Vert x-x_I\Vert \le {\rho }_I,x\in \Omega \right\}$ . For a given point $x\in \Omega$ , define the indicator set ${\tau }_x=\left\{I|\Vert x_I-x\Vert <{\rho }_I,x_I\in X\right\}$ .

This subsection is adopted from the book by Liu and Gu [24]. We consider an unknown scalar function of the field variable $u\left(x\right)$ in the region $\Omega$ . The IMLS approximation of $u\left(x\right)$ at $x$ is defined by $u^h\left(x\right)$ . In order to achieve interpolation property, we adopt a singular weight function to construct the orthogonal basis in reference to Lancaster and Salkauskas [25].

The weight function is defined by

$w\left(x,x_I\right)=\{\begin{array}{ll}{m}_I\left(x\right){\Vert \frac{x-x_I}{{\rho }_I}\Vert }^{-\overline{\alpha }},\hfill & \Vert x-x_I\Vert <{\rho }_I,\hfill \\ 0,\hfill & \text{others},\hfill \end{array}$(3.1)

where the parameter $\overline{\alpha }$ is an even positive integer and $m_I\left(x\right)$ can be chosen to be any weight function used in the MLS approximation [26].

Normalized form for interpolation is constructed by

$\nu \left(x,x_I\right)=\frac{w\left(x,x_I\right)}{{{\displaystyle \sum }}_{J\in {\tau }_x}w\left(x,x_J\right)},$ (3.2)

then the function $\nu \left(x,x_I\right)$ has the interpolation properties. Now we introduce the construction of the shape function. $p\left(x\right)$ is the basis function for the spatial coordinates and $m$ is the number of these basis functions. For $x={\left[{\xi }_1,{\xi }_2\right]}^{\text{T}}$ , we usually choose the following basis functions: $p\left(x\right)=\left[1,{\xi }_1,{\xi }_2,\cdots \right]$ , which are constructed using monomials in Pascal’s triangle to ensure minimal completeness. Let the basis functions $p_1\left(x\right)\equiv 1,p_2\left(x\right),\cdots ,p_m\left(x\right)$ be given. We will generate a new set of basis functions from these. First normalize $p_1\left(x\right)$ such that

${\tilde{p}}_1\left(x,x_I\right)=\frac{1}{\sqrt{{{\displaystyle \sum }}_{I}w\left(x,x_I\right)}}.$(3.3)

We can then generate new basis functions orthogonal to ${\tilde{p}}_1$ :

${\tilde{p}}_i\left(x,x_I\right)=p_i\left(x_I\right)-\mathcal{S}{p}_i\left(x\right),i=2,\cdots ,m,$(3.4)

where $\mathcal{S}$ is a linear operator defined by

$\mathcal{S}{p}_i\left(x\right)={\displaystyle \sum }_{I\in {\tau }_x}\nu \left(x,x_I\right)p_i\left(x_I\right),i=2,\cdots ,m.$ (3.5)

In order to obtain an expression for the approximation function $u^h\left(x\right)$ that satisfies the interpolation property, Lancaster and Sarkauskas [25] defined a local approximation:

$u^h\left(x,x_I\right)={\tilde{p}}_1\left(x,x_I\right)a_1\left(x\right)+{\displaystyle \sum }_{i=2}^m{\tilde{p}}_i\left(x,x_I\right)a_i\left(x\right),$(3.6)

where $a_i\left(x\right) \left(i=1,\cdots ,m\right)$ are the coefficients and can be solved later. For a given $x$ , the difference between the locally approximated function $u^h\left(x,x_I\right)$ and the function $u\left(x_I\right)$ is minimized by weighted least squares. The weighted discrete $L^2$ -norm is defined by

$J={\displaystyle \sum }_{I\in {\tau }_x}w\left(x,x_I\right){\left[u^h\left(x,x_I\right)-u_I\right]}^2,$ (3.7)

where $w\left(x,x_I\right)$ (as shown in Equation (3.1)) is the weight function, $x_I$ ($I\in {\tau }_x$ ) are the points in the support domain of $x$ and $u_I=u\left(x_I\right)$ .

Now define the inner product:

${\left(f,g\right)}_x={\displaystyle \sum }_{I\in {\tau }_x}w\left(x,x_I\right)f\left(x_I\right)g\left(x_I\right),\forall f,g\in C^0\left(\Omega \right),$ (3.8)

where the subscript $x$ denotes a point in $\text{Ω}$ . Then, the corresponding norm at point $x$ is defined by

${\Vert f\Vert }_x={\left[{\displaystyle \sum }_{I\in {\tau }_x}w\left(x,x_I\right)f^2\left(x_I\right)\right]}^{1/2}.$ (3.9)

By minimizing the weighted discrete $L^2$ -norm in Equation (3.7), we have

${\left(u(\cdot )-u^h\left(x,\cdot \right),{\tilde{p}}_1\right)}_x=0,$ (3.10)

and

${\left(u(\cdot )-u^h\left(x,\cdot \right),{\tilde{p}}_i\right)}_x=0,i=2,\cdots ,m.$(3.11)

By orthogonality, Equations (3.10) and (3.11) can be rewritten as

$a_1\left(x\right)={\left(u,{\tilde{p}}_1\right)}_x,$ (3.12)

and

$a_1\left(x\right){\left({\tilde{p}}_1,{\tilde{p}}_j\right)}_x+{\displaystyle \sum }_{i=2}^m{a}_i\left(x\right){\left({\tilde{p}}_i,{\tilde{p}}_j\right)}_x={\left(u,{\tilde{p}}_j\right)}_x,j=2,\cdots ,m.$(3.13)

According to Equations (3.3), (3.5), (3.12) and the definition of the inner product, we have

${\tilde{p}}_1\left(x,x_I\right)a_1\left(x\right)=\frac{1}{{\left[{{\displaystyle \sum }}_{I\in {\tau }_x}w\left(x,x_I\right)\right]}^{1/2}}{\left(u,{\tilde{p}}_1\right)}_x={\displaystyle \sum }_{I\in {\tau }_x}\nu \left(x,x_I\right)u_I=\mathcal{S}u.$(3.14)

Then Equation (3.13) can be simplified as

${\displaystyle \sum }_{i=2}^m{a}_i\left(x\right){\left({\tilde{p}}_i,{\tilde{p}}_j\right)}_x={\left(u-\mathcal{S}u,{\tilde{p}}_j\right)}_x,j=2,\cdots ,m.$(3.15)

In reference [25], the unknown parameter $a_i\left(x\right)$ ($i=2,3,\cdots ,m$ ) is solved by Equation (3.15). In fact, Equation (3.1) can be made even simpler by using the following lemma.

Lemma 5. [26] If the weight function defined by (3.1) is used, there exist

${\left(\mathcal{S}u,{\tilde{p}}_i\right)}_x=0,i=2,\cdots ,m.$

According to Lemma 5, Equation (3.15) can be reduced to

${\displaystyle \sum }_{i=2}^m{a}_i\left(x\right){\left({\tilde{p}}_i,{\tilde{p}}_j\right)}_x={\left(u,{\tilde{p}}_j\right)}_x,j=2,\cdots ,m.$(3.16)

Equation (3.16) is simpler to the corresponding expression in [25] and can be rewritten as

$A\left(x\right)a\left(x\right)=F_W\left(x\right)u,$ (3.17)

where

$u^{\text{T}}=\left(u_1,u_2,\cdots ,u_M\right),$ (3.18)

$a^{\text{T}}\left(x\right)=\left(a_2\left(x\right),a_3\left(x\right),\cdots ,a_m\left(x\right)\right),$ (3.19)

$A\left(x\right)=F_W\left(x\right)F^{\text{T}}\left(x\right),$ (3.20)

$F\left(x\right)=\left[\begin{array}{cccc}{\tilde{p}}_2\left(x,x_1\right)& {\tilde{p}}_2\left(x,x_2\right)& \cdots & {\tilde{p}}_2\left(x,x_M\right)\\ {\tilde{p}}_3\left(x,x_1\right)& {\tilde{p}}_3\left(x,x_2\right)& \cdots & {\tilde{p}}_3\left(x,x_M\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{p}}_m\left(x,x_1\right)& {\tilde{p}}_m\left(x,x_2\right)& \cdots & {\tilde{p}}_m\left(x,x_M\right)\end{array}\right],$ (3.21)

and $F_W\left(x\right)={\left({\varpi }_{kj}\left(x\right)\right)}_{m\times M}$ is a $m\times M$ matrix with

${\varpi }_{kj}\left(x\right)=\{\begin{array}{ll}w\left(x,x_j\right){\tilde{p}}_k\left(x,x_j\right),\hfill & x\ne x_j,\hfill \\ {\displaystyle {\sum }_{I\in {\tau }_x,I\ne j}w}\left(x,x_I\right)\left[{\tilde{p}}_k\left(x_j\right)-{\tilde{p}}_k\left(x_I\right)\right],\hfill & x=x_j.\hfill \end{array}$

Here $M$ is the number of nodes in the support domain.

Then from Equation (3.17) we have

$a\left(x\right)=A^{-1}\left(x\right)F_W\left(x\right)u.$ (3.22)

In turn, bringing in Equation (3.6) yields the local approximation function as

$u^h\left(x,x_I\right)=\mathcal{S}u+{\displaystyle \sum }_{i=2}^m{a}_i\left(x\right){\tilde{p}}_i\left(x,x_I\right).$(3.23)

Then we denote $g_i$ by ${\tilde{p}}_i$ and this yields the global interpolating approximation function for $u\left(x\right)$ :

$u^h\left(x\right)=\mathcal{S}u+{\displaystyle \sum }_{i=2}^m{a}_i\left(x\right)g_i\left(x\right)\equiv \Phi \left(x\right)u={\displaystyle \sum }_{i=1}^M{\Phi }_i\left(x\right)u\left(x_i\right),$ (3.24)

where $\Phi \left(x\right)=\left({\Phi }_1\left(x\right),{\Phi }_2\left(x\right),\cdots ,{\Phi }_M\left(x\right)\right)$ is the shape function matrix. Its expression is

$\Phi \left(x\right)={\nu }^{\text{T}}+p^{\text{T}}\left(x\right)A^{-1}\left(x\right)F_W\left(x\right),$ (3.25)

where the constituent vectors are defined as follows:

${\nu }^{\text{T}}=\left(\nu \left(x,x_1\right),\cdots ,\nu \left(x,x_M\right)\right),$ (3.26)

$p^{\text{T}}=\left(g_2\left(x\right),\cdots ,g_m\left(x\right)\right),$ (3.27)

$g_i=p_i\left(x\right)-\mathcal{S}{p}_i\left(x\right).$ (3.28)

The shape functions:

${\Phi }_I\left(x\right)=\nu \left(x,x_I\right)+{\displaystyle \sum }_{i=2}^m{\tilde{p}}_i\left(x\right){\left(A^{-1}\left(x\right)F_W\left(x\right)\right)}_{iI},I=1,\cdots ,M,$(3.29)

satisfying ${\Phi }_I\left(x_J\right)={\delta }_{IJ}$ .

Furthermore, assuming that $u\in H^{m+1}\left(\Omega \right)$ , the properties of approximate solutions can be obtained [26]:

${\Vert D^{\mu }u-D^{\mu }{u}^h\Vert }_{L^{{}_2}\left(\Omega \right)}\le C{\delta }^{m+1-\left|\mu \right|}{\Vert u\Vert }_{m+1},\mu \le 2,$ (3.30)

where $\delta$ represents the node distance parameter in $\Omega$ , which will be applied extensively in the following sections.

3.2. Details of Constructing the IMLS Shape Function

This subsection employs the theoretical foundation established in the previous subsection to construct the shape functions, taking a 2-dimension (2D) fourth-order equation as a representative example. The equation is introduced as follows:

${\Delta }^{2}u\left({\xi }_1,{\xi }_2\right)-\Delta u\left({\xi }_1,{\xi }_2\right)=f\left({\xi }_1,{\xi }_2\right),\left({\xi }_1,{\xi }_2\right)\in \Omega =\left[0,1\right]\times \left[0,1\right],$

with boundary conditions: $u=0$ and $\text{Δ}u=0$ . Generate 25 uniformly distributed nodes within $\Omega =\left[0,1\right]\times \left[0,1\right]$ :

• Node coordinates: $\left(i/4,j/4\right)$ where $i,j=0,1,2,3,4$ .

• Node numbering: $I=j\times 5+i+1$ (25 nodes in total, numbered 1 - 25).

• Minimum node distance: $h=0.25$ .

• Radius of the influence domain: $\rho =1.5h$ (ensuring the influence domain of each calculation point contains approximately 9 nodes).

First, we construct the cubic polynomial basis functions. The complete set of 2D original cubic polynomial basis functions (including all terms of degree ≤ 3) is:

$\begin{array}{l}{p}_1\left({\xi }_1,{\xi }_2\right)=1,p_2\left({\xi }_1,{\xi }_2\right)={\xi }_1,p_3\left({\xi }_1,{\xi }_2\right)={\xi }_2,p_4\left({\xi }_1,{\xi }_2\right)={\xi }_1^2,\\ p_5\left({\xi }_1,{\xi }_2\right)={\xi }_1{\xi }_2,p_6\left({\xi }_1,{\xi }_2\right)={\xi }_2^2,p_7\left({\xi }_1,{\xi }_2\right)={\xi }_1^3,\\ p_8\left({\xi }_1,{\xi }_2\right)={\xi }_1^2{\xi }_2,p_9\left({\xi }_1,{\xi }_2\right)={\xi }_1{\xi }_2^2,p_{10}\left({\xi }_1,{\xi }_2\right)={\xi }_2^3.\end{array}$

The number of basis functions is $m=10$ . The objective is to determine the method for approximating second-order derivatives in fourth-order equations.

Second, orthogonalization is applied to the basis functions. For the sake of illustration, we will consider the calculation point $\left(0.25,0.25\right)$ .

Step 1: the nodes in the influence domain of the calculation point $x_q=\left(0.25,0.25\right)$ are:

$\begin{array}{l}\left(0,0\right),\left(0.25,0\right),\left(0.5,0\right),\left(0,0.25\right),\left(0.25,0.25\right),\\ \left(0.5,0.25\right),\left(0,0.5\right),\left(0.25,0.5\right),\left(0.5,0.5\right).\end{array}$

There are 9 nodes in total (numbered: 1, 2, 3, 6, 7, 8, 11, 12, 13).

Step 2: an improved weight function is adopted:

$w_I\left(x\right)=w\left(x,x_I\right)={\left(\Vert \frac{x-x_I}{\rho }\Vert \right)}^{-2},\Vert x-x_I\Vert \le \rho .$

In order to avoid numerical singularities, when $x=x_I\left(w=\infty \right)$ we take $w={10}^4$ (a sufficiently large value to replace infinity). Then the calculated weight values for each node are as follows (examples): central node $x_5=\left(0.25,0.25\right)$ : $w_5=+\infty$ (taking 10⁴), adjacent node $x_2=\left(0.25,0\right)$ : distance = 0.25, $w_2={\left(0.25/\rho \right)}^{-2}=2.250$ , and diagonal node $x_1=\left(0,0\right)$ : distance ≈ 0.354, $w_1=1.125$ . Normalization weight function:

${\nu }_I=\frac{w_I}{{{\displaystyle \sum }}_{J=1}^9{w}_J}.$

Step 3: orthogonalization calculation (Gram-Schmidt method):

1) Zero-order Orthogonal Basis ${\tilde{p}}_1$ :

${\tilde{p}}_1=\frac{p_1}{\sqrt{{{\displaystyle \sum }}_I{w}_I{p}_1{\left(x_I\right)}^2}}\approx \frac{1}{\sqrt{{10}^4}}.$

2) First-order Orthogonal Basis ${\tilde{p}}_2$ (Eliminating Correlation with ${\tilde{p}}_1$ ):

${\alpha }_1=\frac{{{\displaystyle \sum }}_I{w}_I{p}_2\left(x_I\right){\tilde{p}}_1\left(x_I\right)}{{{\displaystyle \sum }}_I{w}_I{\tilde{p}}_1{\left(x_I\right)}^2},$

${\tilde{p}}_2=p_2-{\alpha }_1{\tilde{p}}_1={\xi }_1-{\alpha }_1{\tilde{p}}_1.$

3) Subsequent Orthogonal Basis (Recursive Formula): for $k=3,4,\cdots ,10$ :

${\tilde{p}}_k=p_k-{\displaystyle \sum }_{j=1}^{k-1}\frac{{\left(p_k,{\tilde{p}}_j\right)}_x}{{\left({\tilde{p}}_j,{\tilde{p}}_j\right)}_x}{\tilde{p}}_j,$

where inner product is defined by (0.23). Then, we can construct subsequent orthogonal basis functions according to this procedure.

Third, assemble the shape functions. For the calculation point $x_q=\left(0.25,0.25\right)$ , the shape functions ${\Phi }_I$ are obtained through the following steps:

Step 1: construct the basis function matrix ${\left(F\right)}_{10\times 9}$ : each row corresponds to an orthogonal basis, each column corresponds to a node within the influence domain, and the values are ${\tilde{p}}_k\left(x_I\right)$ $\left(k=1,\cdots ,m\right)$ .

Step 2: construct the weight matrix ${\left(W\right)}_{9\times 9}$ (diagonal matrix): The diagonal elements are the weight function values $w_I$ , and define $F_W=FW$ .

Step 3: solve the coefficient matrix:

${\left(A\right)}_{10\times 10}=F_W{F}^{\text{T}},{\left(\lambda \right)}_{10\times 9}=A^{-1}PW,a=A^{-1}{F}_{W}U,$

where $U$ represents the function values at the 25 nodes.

Step 4: calculate the shape functions:

${\Phi }_I\left(x_q\right)={\nu }^{\text{T}}+p^{\text{T}}\left(x_q\right)A^{-1}{F}_W,$

where ${\nu }^{\text{T}}=\left[{\nu }_1,\cdots ,{\nu }_9\right]$ and $p^{\text{T}}\left(x_q\right)=\left[{\tilde{p}}_1\left(x_q\right),{\tilde{p}}_2\left(x_q\right),\cdots ,{\tilde{p}}_{10}\left(x_q\right)\right]$.

Finally, 9 non-zero shape functions are obtained, satisfying:

• Interpolation property: ${\Phi }_I\left(x_J\right)={\delta }_{IJ}$ (exact interpolation at nodes).

• Continuity of second-order derivatives: both $\frac{{\partial }^2{\Phi }_I}{\partial {\xi }_1^2}$ and $\frac{{\partial }^2{\Phi }_I}{\partial {\xi }_2^2}$ are linear functions.

Now we can construct the discrete scheme for fourth-order equation. Using the weak form ${\displaystyle {\int }_{\Omega }\Delta u\cdot \Delta v\text{d}{\xi }_1\text{d}{\xi }_2}+{\displaystyle {\int }_{\Omega }\nabla u\cdot \nabla v\text{d}{\xi }_1\text{d}{\xi }_2}={\displaystyle {\int }_{\Omega }fv\text{d}{\xi }_1\text{d}{\xi }_2}$ , the equation is discretized as:

$\left(K_1+K_2\right)U=f.$

where:

• Stiffness matrix: $K_1{}_{\left(IJ\right)}={\displaystyle {\int }_{\Omega }\Delta {\Phi }_I\cdot \Delta {\Phi }_J\text{d}{\xi }_1\text{d}{\xi }_2}$ , $K_2{}_{\left(IJ\right)}={\displaystyle {\int }_{\Omega }\nabla {\Phi }_I\cdot \nabla {\Phi }_J\text{d}{\xi }_1\text{d}{\xi }_2}$ .

• Load vector: $f_I={\displaystyle {\int }_{\Omega }f{\Phi }_I\text{d}{\xi }_1\text{d}{\xi }_2}$ .

• Unknown vector: $U$ represents the function values at the 25 nodes.

It is noteworthy that the background integration grid was employed for the calculation of the inner product during the process of matrix assembly.

3.3. Fully Discrete Approximation Schemes for the Optimal Control Problem

The IMLS approximation technique is incorporated into the (OCP-OPT) model derived from the Galerkin weak formulation. By employing the backward-Euler method for temporal discretization, the IEFG method is ultimately developed. The shape function space is $V_J=\text{span}\left\{{\Phi }_1\left(x\right),\cdots ,{\Phi }_M\left(x\right)\right\}$ .

Define the finite dimensional subspace of the space of control variables $K$ :

$K_h=L^2\left(0,T;K\cap V_J\right).$

Define the finite dimensional subspace of the state space $V$ :

$V_h=L^2\left(0,T;V_J\right).$

Define the time discrete target functional for $i=1,\cdots ,N$ :

$J_{hk}\left(v_h^i\right)=\frac{1}{2}{\displaystyle \sum }_{i=1}^N{k}_i\left({\Vert Y_h^i-y_d^i\Vert }_{0,\Omega }^2+\alpha {\Vert v_h^i\Vert }_{0,{\Omega }_U}^2\right),\forall \left(Y_h^i,v_h^i\right)\in V_h\times K_h,$

satisfying

$J_{hk}\left(U_h^i\right)=\underset{v_h^i\in K_h}{\text{min}}{J}_{hk}\left(v_h^i\right).$

It is evident that the control variable $u$ manifests as a control term on the right-hand side of the equation. Consequently, the accuracy requirements for $u$ are less than those for the state variable $y$ and the co-state variable $p$ . This paper proposes a hybrid discrete scheme, which is a method for approximating continuous functions with a combination of different types of functions. In this scheme, the state variables and co-state variables are discretized using IMLS shape functions. These functions satisfy the high-order continuity requirements of fourth-order problems. In the context of the control variable $u$ , the discrete scheme adopts a piecewise constant space, which is derived from background integration grids. This discrete scheme aligns with the principles of the finite element method, ensuring computational accuracy while enhancing computational efficiency and facilitating theoretical analysis.

Now, in order to obtain the fully discrete approximation scheme, let

$\{\begin{array}{l}{U}_h\left(x\right)={\displaystyle {\sum }_{i=1}^{M}u}\left(x_i\right){\psi }_i\left(x\right),\\ Y_h\left(x\right)={\displaystyle {\sum }_{i=1}^{M}y}\left(x_i\right){\varphi }_i\left(x\right),\\ P_h\left(x\right)={\displaystyle {\sum }_{i=1}^{M}p}\left(x_i\right){\varphi }_i\left(x\right),\end{array}$ (3.31)

where ${\varphi }_i$ is the basis function of the IMLS approximation, and ${\psi }_i$ is the basis function of the piecewise constant space derived from background integration grids of meshless method.

This gives the fully discrete approximation scheme: look for $\left(Y_h^i,P_h^{i-1},U_h^i\right)\in \left(V_h\times V_h\times K_h\right)$ $\left(i=1,\cdots ,N\right)$ satisfying

$\{\begin{array}{l}\left(\frac{Y_h^i-Y_h^{i-1}}{k_i},w_h\right)+\left(\Delta Y_h^i,\Delta w_h\right)+a\left(Y_h^i,w_h\right)=\left(f^i+B{U}_h^i,w_h\right),\hfill \\ Y_h^0\left(x\right)=y_0^h\left(x\right),w_h\in V_h,\hfill \end{array}$ (3.32)

$\{\begin{array}{l}\left(\frac{P_h^{i-1}-P_h^i}{k_i},q_h\right)+\left(\Delta P_h^{i-1},\Delta q_h\right)+a\left(P_h^{i-1},q_h\right)=\left(Y_h^i-y_d^i,q_h\right),\hfill \\ P_h^N\left(x\right)=0,q_h\in V_h,\hfill \end{array}$ (3.33)

${\left(\alpha U_h^i+B^{\text{*}}{P}_h^{i-1},v_h-U_h^i\right)}_U\ge 0,\forall v_h\in K_h.$ (3.34)

4. Error Estimates

This section presents a complete derivation of the a prior error estimates for the IEFG approximation of the optimal control problem previously referenced. Prior to embarking on the proof, it is important to introduce some key projections and their associated properties, ensuring a streamlined and effective subsequent argumentation.

4.1. Key Projections

For the sake of subsequent argumentation, we provide the following definitions and related properties of projections at this subsection. First, we introduce the concept of a bilinear form. As mentioned above, $\Omega$ is a bounded Lipschitz domain. Consider the bilinear form defined on the space $V$ :

$A\left(u,v\right)=\left(\Delta u,\Delta v\right)+a\left(u,v\right),\forall u,v\in V.$

We now prove the coercivity and boundedness of this bilinear form.

Theorem 6 (Boundedness) There exists a constant $M>0$ such that

$\left|A\left(u,v\right)\right|\le M{\Vert u\Vert }_{H^2\left(\Omega \right)}{\Vert v\Vert }_{H^2\left(\Omega \right)},\forall u,v\in H_0^2\left(\Omega \right).$

Proof. By the definition of $A\left(u,v\right)$ and the triangle inequality:

$\left|A\left(u,v\right)\right|=\left|\left(\text{Δ}u,\text{Δ}v\right)+a\left(u,v\right)\right|\le \left|\left(\text{Δ}u,\text{Δ}v\right)\right|+\left|a\left(u,v\right)\right|.$

We estimate each term separately. For the first term, by the Cauchy-Schwarz inequality

$\left|\left(\Delta u,\Delta v\right)\right|\le {\Vert \Delta u\Vert }_{L^2\left(\Omega \right)}{\Vert \Delta v\Vert }_{L^2\left(\Omega \right)}\le {\Vert u\Vert }_{H^2\left(\Omega \right)}{\Vert v\Vert }_{H^2\left(\Omega \right)}.$

For the second term, using the boundedness of $a\left(u,v\right)$ and the fact that ${\Vert \cdot \Vert }_{H^1\left(\Omega \right)}\le {\Vert \cdot \Vert }_{H^2\left(\Omega \right)}$ , we can derive