# Solution of the Optimal Control Problem for the Cahn-Hilliard Equation Using
Finite Difference Approximation

Gobinda Garai & Bankim C. Mandal  Gobinda Garai
School of Basic Sciences, Indian Institute of Technology Bhubaneswar, India
<EMAIL_ADDRESS>Bankim C. Mandal
School of Basic Sciences, Indian Institute of Technology Bhubaneswar, India
<EMAIL_ADDRESS>
###### Abstract.

This paper is concerned with the designing, analyzing and implementing linear
and nonlinear discretization scheme for the distributed optimal control
problem (OCP) with the Cahn-Hilliard (CH) equation as constrained. We propose
three difference schemes to approximate and investigate the solution behaviour
of the OCP for the CH equation. We present the convergence analysis of the
proposed discretization. We verify our findings by presenting numerical
experiments.

###### Key words and phrases:

Optimal control problem, Constrained Optimization, Convergence analysis, Cahn-
Hilliard equation.

######  Mathematics Subject Classification:

65M06, 65M12, 49M41

Submitted.

##  1\. Introduction

In this work, we consider the following distributed optimal control problem
(OCP)

$\min\limits_{u(x,t)\in
L^{2}(0,T;L^{2}(\Omega))}J(y,u)=\frac{1}{2}\int_{0}^{T}\int_{\Omega}\left(y(x,t)-\widehat{y}(x,t)\right)^{2}dxdt+\frac{\lambda}{2}\int_{0}^{T}\int_{\Omega}u(x,t)^{2}dxdt,$
(1.1)

subject to the Cahn-Hilliard (CH) equation

$\begin{cases}\frac{\partial y}{\partial t}=\Delta
f(y)-\epsilon^{2}\Delta^{2}y+u,&(x,t)\in\Omega\times(0,T],\\\
\partial_{\nu}y=0=\partial_{\nu}(\Delta
y),&(x,t)\in\partial\Omega\times(0,T],\\\ u(x,0)=u_{0},\end{cases}$ (1.2)

where $\Omega\subset\mathbb{R}^{d}(d=1,2)$ is a bounded Lipschitz domain,
$\nu$ is the unit outward normal, $u\in L^{2}(0,T;L^{2}(\Omega))$ is the
control variable, $\widehat{y}\in L^{2}(0,T;L^{2}(\Omega))$ is the target
state, $\lambda>0$ is a regularization parameter, $0<\epsilon\ll 1$ and
$f(y)=y^{3}-y$. The purpose of the optimization is to determine the control
function $u(x,t)$ in such a way that the state $y(x,t)$ gets as closely as
possible to a given desired state $\widehat{y}(x,t)$. Here we follow optimize-
then-discretize approach to discretize the OCP (1.1) - (1.2).

The existence and uniqueness of OCP (1.1) - (1.2) follows from the standard
variational arguments [14, 18], which can be described by its first order
optimality system, see [Troltzsch]. The first order optimality system yields
the state equation (1.2), the adjoint equation

$\begin{cases}-\frac{\partial p}{\partial t}=f^{\prime}(y)\Delta
p-\epsilon^{2}\Delta^{2}p+\widehat{y}-y,&(x,t)\in\Omega\times(0,T],\\\
\partial_{\nu}p=0=\partial_{\nu}(\Delta
p),&(x,t)\in\partial\Omega\times(0,T],\\\ p(x,T)=0,\end{cases}$ (1.3)

and the optimality condition

$p(x,t)=\lambda u(x,t),(x,t)\in\Omega\times(0,T].$ (1.4)

The CH equation, represented by Equation (1) without the control term $u$, is
a mathematical model commonly used to describe the evolution of a binary
melted alloy below its critical temperature as referenced in [1, 2]. Numerous
research studies have been devoted to developing numerical schemes for
approximating the solution of the CH equation. These schemes often employ
finite difference or finite element methods, with either Dirichlet boundary
conditions [5, 10] or Neumann boundary conditions [8, 11, 16, 9, 7, 17, 4, 3],
among others. For further insights into the numerical approaches used for
solving the CH equation can be found in [15].

The solution of the CH equation encompasses two separate dynamics: phase
separation, characterized by rapid changes over time, and phase coarsening,
which occurs at a slower pace. During the initial stages, the formation of
fine-scale phase regions takes place, with interfaces of width $\epsilon$
separating them. On the other hand, during phase coarsening, the solution
tends to converge towards an equilibrium state. The OCP associated with the CH
equation becomes crucial in situations where there is a desire to exert
influence over the phase separation or coarsening behaviour, or to attain
specific concentration profiles. Within this context, the OCP plays an
important role by providing a framework to actively control and shape the
dynamics of phase separation or coarsening in order to achieve desired
outcomes. The following assumptions are made regarding the nonlinear terms:

$\max_{y}|f^{\prime}(y)|\leq
M,\;\max_{y}|\widetilde{f}^{\prime}(y)|\leq\widetilde{M}$ (1.5)

where $M,\widetilde{M}$ are non-negative constant and
$\widetilde{f}(y)=y^{3}$.

We now talk about a few studies that relate OCP to the CH equation. The
authors of [19, 13] took into account the OCP for the CH equation and
presented the existence results. The studies in [20, 6] are concerned with OCP
for the viscous CH equation and their existence results. You may see a study
for OCP with convective CH as constraints in [21]. The OCP for the CH equation
with state constraint has been studied in [22]. However, the above stated body
of works do not have any numerical scheme to solve the OCP. This is, as far as
we are aware, the first attempt at constructing a numerical framework for the
OCP of the CH equation. In order to discretize continuous spatial and temporal
variables, we employ the finite difference approach.

Rest of the paper is organised as follows. In Section 2 we present the
numerical schemes to approximate the solution. In Section 3 we show the
convergence results of the proposed schemes. Finally we present numerical
results in Section 4.

##  2\. Notation and Difference Schemes

We propose the finite difference scheme for (1.2)-(1.3) in one spatial
dimension, while an extension to higher spatial dimension is straight forward.
To discretize the problem (1.2)-(1.3) in the spatial domain $\Omega=(a,b)$, we
partitioned the domain uniformly by taking mesh size $h=(b-a)/(N-1)$, where
$N$ is the number of spatial grid point. Then the discrete domain is
$\Omega_{h}=\\{(x_{i-1},x_{i}):x_{i}=a+(i-1)h,i=1,2,...,N\\}.$ Let
$\delta_{t}=T/N_{t}$ be the uniform time step corresponding to discrete
temporal variable $t_{n}=n\delta_{t}$ for $n=0,1,...,N_{t}$, where $N_{t}$ is
any positive integer. For any function $z$ define on $\Omega_{h}$, we denote
$z_{i}=z(x_{i})$ and define the following difference operators

$\nabla_{h}^{+}z_{i}=\frac{z_{i+1}-z_{i}}{h},\;\nabla_{h}^{-}z_{i}=\frac{z_{i}-z_{i-1}}{h},\;\Delta_{h}=\nabla_{h}^{+}\nabla_{h}^{-}=\nabla_{h}^{-}\nabla_{h}^{+},\;\Delta_{h}^{2}=\Delta_{h}\Delta_{h}.$

Now we define the discrete $L^{2}$-inner product as
$(y,z)_{h}=\sum\limits_{i=1}^{N}y_{i}z_{i}$ for any function $y,z$, defined on
$\Omega_{h}$. Then the corresponding discrete $L^{2}$-norm is given by
$\|z\|_{h}=\sqrt{(z,z)_{h}}$. We also define discrete $H^{1}$-seminorm,
$H^{2}$-seminorm and maximum-norm for any grid function $z$ in $\Omega_{h}$ as

$|z|_{1,h}^{2}=h\sum\limits_{i=1}^{N}(\nabla_{h}^{-}z_{i})^{2},\;|z|_{2,h}^{2}=h\sum\limits_{i=1}^{N}(\Delta_{h}z_{i})^{2},\;\|z\|_{\infty,h}=\sup\limits_{i}|z_{i}|.$

Let $Y_{i}^{n}=Y(x_{i},t_{n}),P_{i}^{n}=P(x_{i},t_{n})$ and
$U_{i}^{n}=U(x_{i},t_{n})$ for $i=1,2,...,N$ and $n=0,1,...,N_{t}$. Then to
get the approximate solution of (1.2)-(1.3) we discretize the state equation
(1.2) as

$\frac{Y_{i}^{n+1}-Y_{i}^{n}}{\delta_{t}}=\Delta_{h}f(Y_{i}^{n+1})-\epsilon^{2}\Delta_{h}^{2}Y_{i}^{n+1}+\frac{1}{\lambda}P_{i}^{n},$
(2.1)

and the adjoint equation (1.3) as

$-\frac{P_{i}^{n+1}-P_{i}^{n}}{\delta_{t}}=f^{\prime}(Y_{i}^{n})\Delta_{h}P_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}P_{i}^{n}+\widehat{Y}_{i}^{n+1}-Y_{i}^{n+1}.$
(2.2)

Observe that we made a substitution of control variable $U_{i}^{n}$ by
$\frac{1}{\lambda}P_{i}^{n}$ in (2.1) using the relation in (1.4). The
described scheme (2.1) - (2.2) is undoubtedly nonlinear, with the nonlinearity
originating from (2.1). We call the scheme (2.1) - (2.2) as $\textbf{S}1$.
Next we display two possible linear approximation to the state equation (1.2).
First we consider the following linear scheme for the state equation

$\frac{Y_{i}^{n+1}-Y_{i}^{n}}{\delta_{t}}=\Delta_{h}((Y_{i}^{n})^{2}Y_{i}^{n+1})-\Delta_{h}Y_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}Y_{i}^{n+1}+\frac{1}{\lambda}P_{i}^{n}.$
(2.3)

We call the state approximation scheme (2.3), along with adjoint approximation
scheme (2.2) as $\textbf{S}2$. Another way to discretize the state equation is

$\frac{Y_{i}^{n+1}-Y_{i}^{n}}{\delta_{t}}=\Delta_{h}(Y_{i}^{n})^{3}-3\Delta_{h}Y_{i}^{n}+2\Delta_{h}Y_{i}^{n+1}-\epsilon^{2}\Delta_{h}^{2}Y_{i}^{n+1}+\frac{1}{\lambda}P_{i}^{n}.$
(2.4)

We call the state approximation scheme (2.4), along with adjoint approximation
scheme (2.2) as $\textbf{S}3$.

##  3\. Convergence of the Schemes

We examine the convergence of the finite difference scheme (2.1)-(2.2) in this
section. Discrete Gronwall’s inequality, which will be used extensively in our
convergence estimations, was previously discussed in [12].

######  Lemma 3.1 (Discrete Gronwall inequality).

Let $\tau,Q$ and $a_{n},b_{n},c_{n},d_{n}$ be non-negative numbers for
integers $n\geq 0$ such that

$a_{n}+\tau\sum\limits_{k=0}^{n}b_{k}\leq\tau\sum\limits_{k=0}^{n}d_{k}a_{k}+\tau\sum\limits_{k=0}^{n}c_{k}+Q.$

Suppose that $\tau d_{k}<1,\forall k$ and set $\sigma_{k}=1/(1-\tau d_{k})$.
Then the following holds for $n\geq 0$

$a_{n}+\tau\sum\limits_{k=0}^{n}b_{k}\leq\exp\left(\tau\sum\limits_{k=1}^{n}\sigma_{k}d_{k}\right)\left[\tau\sum\limits_{k=0}^{n}c_{k}+Q\right].$

######  Theorem 3.2 (Convergence of the scheme $\textbf{S}1$).

Let $y(x,t)$ and $p(x,t)$ be sufficiently smooth functions. For $\delta_{t}$
sufficiently small, the finite difference scheme (2.1)-(2.2) is first order in
time and second order in space convergent, i.e.,

$\max\limits_{n}\left\\{\parallel y^{n}-Y^{n}\parallel_{\infty,h}+\parallel
p^{n}-P^{n}\parallel_{\infty,h}\right\\}\leq C(\delta_{t}+h^{2}).$

###### Proof.

Using Taylor expansion we observe that the exact solution
$y_{i}^{n}=y(x_{i},t_{n})$ and $p_{i}^{n}=p(x_{i},t_{n})$ satisfy the
following equations

$\frac{y_{i}^{n+1}-y_{i}^{n}}{\delta_{t}}=\Delta_{h}f(y_{i}^{n+1})-\epsilon^{2}\Delta_{h}^{2}y_{i}^{n+1}+\frac{1}{\lambda}p_{i}^{n}+F_{i}^{n},$
(3.1)
$-\frac{p_{i}^{n+1}-p_{i}^{n}}{\delta_{t}}=f^{\prime}(y_{i}^{n})\Delta_{h}p_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}p_{i}^{n}+\widehat{y}_{i}^{n+1}-y_{i}^{n+1}+G_{i}^{n},$
(3.2)

where $F_{i}^{n}$ and $G_{i}^{n}$ denote the truncation error, which satisfy
the following for some positive constants $c_{1},c_{2}$

$\max\limits_{i,n}|F_{i}^{n}|\leq
c_{1}(\delta_{t}+h^{2}),\;\;\max\limits_{i,n}|G_{i}^{n}|\leq
c_{2}(\delta_{t}+h^{2}).$ (3.3)

Let us define the error $e_{i}^{n}=y_{i}^{n}-Y_{i}^{n}$ and
$\mathfrak{e}_{i}^{n}=p_{i}^{n}-P_{i}^{n}$. Then taking difference between
(3.1)-(3.2) and (2.1)-(2.2) yields

$\frac{e_{i}^{n+1}-e_{i}^{n}}{\delta_{t}}=\Delta_{h}\left(f(y_{i}^{n+1})-f(Y_{i}^{n+1})\right)-\epsilon^{2}\Delta_{h}^{2}e_{i}^{n+1}+\frac{1}{\lambda}\mathfrak{e}_{i}^{n}+F_{i}^{n},$
(3.4)
$-\frac{\mathfrak{e}_{i}^{n+1}-\mathfrak{e}_{i}^{n}}{\delta_{t}}=f^{\prime}(y_{i}^{n})\Delta_{h}p_{i}^{n}-f^{\prime}(Y_{i}^{n})\Delta_{h}P_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}\mathfrak{e}_{i}^{n}-e_{i}^{n+1}+G_{i}^{n}.$
(3.5)

Taking the inner product of (3.4) and $e^{n+1}$ yields

$\left(\frac{e^{n+1}-e^{n}}{\delta_{t}},e^{n+1}\right)_{h}=\left(\Delta_{h}\left(f(y^{n+1})-f(Y^{n+1})\right),e^{n+1}\right)_{h}-\epsilon^{2}|e^{n+1}|_{2,h}^{2}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},e^{n+1}\right)_{h}+\left(F^{n},e^{n+1}\right)_{h}.$
(3.6)

Using the differentiability of $f$, Cauchy-Schwarz inequality and
$\left(\frac{e^{n+1}-e^{n}}{\delta_{t}},e^{n+1}\right)_{h}\geq\frac{\|e^{n+1}\|_{h}^{2}-\|e^{n}\|_{h}^{2}}{2\delta_{t}}$
on (3.6) we obtain

$\frac{\|e^{n+1}\|_{h}^{2}-\|e^{n}\|_{h}^{2}}{2\delta_{t}}\leq
M\|e^{n+1}\|_{h}|e^{n+1}|_{2,h}-\epsilon^{2}|e^{n+1}|_{2,h}^{2}+\frac{1}{\lambda}\|\mathfrak{e}^{n}\|_{h}\|e^{n+1}\|_{h}+\|F^{n}\|_{h}\|e^{n+1}\|_{h}.$
(3.7)

Using Young’s inequality and (3.3) on (3.7) we get

$\frac{\|e^{n+1}\|_{h}^{2}-\|e^{n}\|_{h}^{2}}{2\delta_{t}}\leq\underbrace{\left(\frac{M^{2}}{4\epsilon^{2}}+\frac{1}{2\lambda}+\frac{1}{2}\right)}_{\alpha}\|e^{n+1}\|_{h}^{2}+\frac{1}{2\lambda}\|\mathfrak{e}^{n}\|_{h}^{2}+\frac{c_{1}}{2}(\delta_{t}+h^{2})^{2}.$
(3.8)

From (3.8) we have

$\left(1-2\delta_{t}\alpha\right)\parallel
e^{n+1}\parallel_{h}^{2}\leq\parallel
e^{n}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+c_{1}\delta_{t}(\delta_{t}+h^{2})^{2}.$
(3.9)

For $\delta_{t}<\frac{1}{2\alpha}$ in (3.9) and taking the sum over $n$ we
have

$\parallel
e^{n}\parallel_{h}^{2}\leq\frac{1}{\left(1-2\delta_{t}\alpha\right)}\sum\limits_{i=1}^{n}\parallel
e^{i}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda\left(1-2\delta_{t}\alpha\right)}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+\frac{c_{1}T}{\left(1-2\delta_{t}\alpha\right)}(\delta_{t}+h^{2})^{2}.$
(3.10)

An application of Gronwall’s inequality for sufficiently small $\delta_{t}$
yields

$\parallel e^{n}\parallel_{h}^{2}\leq
C\delta_{t}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.11)

Taking now the inner product of (3.5) and $\mathfrak{e}^{n}$ yields

$\left(-\frac{\mathfrak{e}^{n+1}-\mathfrak{e}^{n}}{\delta_{t}},\mathfrak{e}^{n}\right)_{h}=\left(f^{\prime}(y^{n})\Delta_{h}p^{n}-f^{\prime}(Y^{n})\Delta_{h}P^{n},\mathfrak{e}^{n}\right)_{h}-\epsilon^{2}|\mathfrak{e}^{n+1}|_{2,h}^{2}-\left(e^{n+1},\mathfrak{e}^{n}\right)_{h}+\left(G^{n},\mathfrak{e}^{n}\right)_{h}.$
(3.12)

Note that

$\left(f^{\prime}(y^{n})\Delta_{h}p^{n}-f^{\prime}(Y^{n})\Delta_{h}P^{n},\mathfrak{e}^{n}\right)_{h}\leq\left(|\left(f^{\prime}(y^{n})+f^{\prime}(Y^{n})\right)|\,|\Delta_{h}\mathfrak{e}^{n}|,\mathfrak{e}^{n}\right)_{h}\leq
2M\|\mathfrak{e}^{n}\|_{h}|\mathfrak{e}^{n+1}|_{2,h}.$ (3.13)

Using (3.13) and
$\left(-\frac{\mathfrak{e}^{n+1}-\mathfrak{e}^{n}}{\delta_{t}},\mathfrak{e}^{n}\right)_{h}\geq\frac{\parallel\mathfrak{e}^{n}\parallel_{h}^{2}-\parallel\mathfrak{e}^{n+1}\parallel_{h}^{2}}{2\delta_{t}}$
on (3.12) we obtain

$\frac{\parallel\mathfrak{e}^{n}\parallel_{h}^{2}-\parallel\mathfrak{e}^{n+1}\parallel_{h}^{2}}{2\delta_{t}}\leq
2M\parallel\mathfrak{e}^{n}\parallel_{h}|\mathfrak{e}^{n+1}|_{2,h}-\epsilon^{2}|\mathfrak{e}^{n+1}|_{2,h}^{2}-\left(e^{n+1},\mathfrak{e}^{n}\right)_{h}+\left(G^{n},\mathfrak{e}^{n}\right)_{h}.$
(3.14)

Using Young’s inequality and (3.3) we get

$\frac{\parallel\mathfrak{e}^{n}\parallel_{h}^{2}-\parallel\mathfrak{e}^{n+1}\parallel_{h}^{2}}{2\delta_{t}}\leq\underbrace{\left(\frac{M^{2}}{\epsilon^{2}}+1\right)}_{\beta}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{1}{2}\parallel{e}^{n+1}\parallel_{h}^{2}+\frac{c_{2}}{2}(\delta_{t}+h^{2})^{2}.$
(3.15)

Thus we obtain

$(1-2\delta_{t}\beta)\parallel\mathfrak{e}^{n}\parallel_{h}^{2}\leq\parallel\mathfrak{e}^{n+1}\parallel_{h}^{2}+\delta_{t}\parallel{e}^{n+1}\parallel_{h}^{2}+c_{2}\delta_{t}(\delta_{t}+h^{2})^{2}.$
(3.16)

For $\delta_{t}<\frac{1}{2\beta}$ in (3.16) and taking the sum over $n$ we
have

$\parallel\mathfrak{e}^{n}\parallel_{h}^{2}\leq\frac{1}{(1-2\delta_{t}\beta)}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+\frac{\delta_{t}}{(1-2\delta_{t}\beta)}\sum\limits_{i=1}^{n}\parallel{e}^{i}\parallel_{h}^{2}+\frac{c_{2}T}{(1-2\delta_{t}\beta)}(\delta_{t}+h^{2})^{2}.$
(3.17)

Applying Gronwall’s inequality on (3.17) for sufficiently small $\delta_{t}$
we obtain

$\parallel\mathfrak{e}^{n}\parallel_{h}^{2}\leq
C\delta_{t}\sum\limits_{i=1}^{n}\parallel
e^{i}\parallel_{h}^{2}+C(\delta_{t}+h^{2})^{2}.$ (3.18)

By adding (3.11) and (3.18) we have

$\parallel
e^{n}\parallel_{h}^{2}+\parallel\mathfrak{e}^{n}\parallel_{h}^{2}\leq
C\delta_{t}\sum\limits_{i=1}^{n}\parallel
e^{i}\parallel_{h}^{2}+C\delta_{t}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.19)

Again using Gronwall’s inequality on (3.19) for sufficiently small
$\delta_{t}$ we get

$\parallel e^{n}\parallel_{h}+\parallel\mathfrak{e}^{n}\parallel_{h}\leq
C(\delta_{t}+h^{2}).$ (3.20)

Next we take the inner product of (3.4) and
$\partial_{t}e^{n}:=\frac{e^{n+1}-e^{n}}{\delta_{t}}$, which yields

$\parallel\partial_{t}e^{n}\parallel^{2}_{h}=\left(\Delta_{h}\left(f(y^{n+1})-f(Y^{n+1})\right),\partial_{t}e^{n}\right)_{h}-\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},\partial_{t}e^{n}\right)_{h}+\left(F^{n},\partial_{t}e^{n}\right)_{h}.$
(3.21)

Using the Lipschitz condition of $f$ and Cauchy-Schwarz inequality on (3.21)
we get

$\parallel\partial_{t}e^{n}\parallel^{2}_{h}+\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}\leq
M\parallel
e^{n+1}\parallel_{h}\parallel\partial_{t}(\Delta_{h}e^{n})\parallel_{h}\\\
+\frac{1}{\lambda}\parallel\mathfrak{e}^{n}\parallel_{h}\parallel\partial_{t}e^{n}\parallel_{h}+\parallel
F^{n}\parallel_{h}\parallel\partial_{t}e^{n}\parallel_{h}.$ (3.22)

Using Young’s inequality and
$\left(\partial_{t}(\Delta_{h}e^{n}),\Delta_{h}e^{n+1}\right)_{h}\geq\frac{|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}}{2\delta_{t}}$
on (3.22) we have the following

$\epsilon^{2}\frac{|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}}{2\delta_{t}}\leq\frac{M^{2}}{\epsilon^{2}}\parallel
e^{n+1}\parallel_{h}^{2}+\frac{\epsilon^{2}}{4}\parallel\partial_{t}(\Delta_{h}e^{n})\parallel_{h}^{2}+\frac{1}{2\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{1}{2}\parallel
F^{n}\parallel_{h}^{2}.$ (3.23)

Using (3.20) and (3.3) on (3.23) we have

$\epsilon^{2}\frac{|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}}{2\delta_{t}}\leq\frac{\epsilon^{2}}{4}\frac{|e^{n+1}|_{2,h}^{2}+|e^{n}|_{2,h}^{2}}{\delta_{t}}+C(\delta_{t}+h^{2})^{2}.$
(3.24)

Thus the equation (3.24) implies

$\frac{\epsilon^{2}}{2}|e^{n+1}|_{2,h}^{2}\leq\frac{3\epsilon^{2}}{2}|e^{n}|_{2,h}^{2}+2\delta_{t}C(\delta_{t}+h^{2})^{2}.$
(3.25)

Taking the sum over $n$ on (3.25) we find that

$\frac{\epsilon^{2}}{2}|e^{n}|_{2,h}^{2}\leq\frac{3\epsilon^{2}}{2}\sum\limits_{i=1}^{n}|e^{i}|_{2,h}^{2}+2TC(\delta_{t}+h^{2})^{2}.$
(3.26)

Applying Gronwall’s inequality on (3.26) we get

$|e^{n}|_{2,h}\leq C(\delta_{t}+h^{2}).$ (3.27)

Also we have
$|e^{n}|_{1,h}^{2}=-(\Delta_{h}e^{n},e^{n})\leq|e^{n}|_{2,h}\parallel
e^{n}\parallel_{h}$. Then using (3.27) and (3.20) we obtain

$|e^{n}|_{1,h}\leq C(\delta_{t}+h^{2}).$ (3.28)

Now we take the inner product of (3.5) and $-\partial_{t}\mathfrak{e}^{n}$,
that yields

$\parallel\partial_{t}\mathfrak{e}^{n}\parallel_{h}^{2}=\left(f^{\prime}(y^{n})\Delta_{h}p^{n}-f^{\prime}(Y^{n})\Delta_{h}P^{n},-\partial_{t}\mathfrak{e}^{n}\right)_{h}-\epsilon^{2}\left(\Delta_{h}\mathfrak{e}^{n+1},-\partial_{t}(\Delta_{h}\mathfrak{e}^{n})\right)_{h}\\\
-\left(e^{n+1},-\partial_{t}\mathfrak{e}^{n}\right)_{h}+\left(G^{n},-\partial_{t}\mathfrak{e}^{n}\right)_{h}.$
(3.29)

Following (3.13) and using Cauchy-Schwarz inequality on (3.29) we have

$\parallel\partial_{t}\mathfrak{e}^{n}\parallel_{h}^{2}+\epsilon^{2}\left(\Delta_{h}\mathfrak{e}^{n+1},-\partial_{t}(\Delta_{h}\mathfrak{e}^{n})\right)_{h}\leq
2M|\mathfrak{e}^{n}|_{2,h}\parallel\partial_{t}\mathfrak{e}^{n}\parallel_{h}+\parallel
e^{n+1}\parallel_{h}\parallel\partial_{t}\mathfrak{e}^{n}\parallel_{h}+\parallel
G^{n}\parallel_{h}\parallel\partial_{t}\mathfrak{e}^{n}\parallel_{h}.$ (3.30)

Using Young’s inequality on (3.30) we get

$\epsilon^{2}\frac{|\mathfrak{e}^{n}|_{2,h}^{2}-|\mathfrak{e}^{n+1}|_{2,h}^{2}}{2\delta_{t}}\leq
4M^{2}|\mathfrak{e}^{n}|_{2,h}^{2}+\frac{1}{2}\parallel
e^{n+1}\parallel_{h}^{2}+\parallel G^{n}\parallel_{h}^{2}.$ (3.31)

Using (3.3) and (3.20) on (3.31) we have

$(\epsilon^{2}-8M^{2}\delta_{t})|\mathfrak{e}^{n}|_{2,h}^{2}\leq\epsilon^{2}|\mathfrak{e}^{n+1}|_{2,h}^{2}+2\delta_{t}C(\delta_{t}+h^{2})^{2}.$
(3.32)

Letting $\delta_{t}<\frac{\epsilon^{2}}{8M^{2}}$ in (3.32) and summing over
$n$ yields

$|\mathfrak{e}^{n}|_{2,h}^{2}\leq\frac{\epsilon^{2}}{(\epsilon^{2}-8M^{2}\delta_{t})}\sum\limits_{i=0}^{n-1}|\mathfrak{e}^{i}|_{2,h}^{2}+2TC(\delta_{t}+h^{2})^{2}.$
(3.33)

Employing Gronwall’s inequality in (3.33) we obtain

$|\mathfrak{e}^{n}|_{2,h}\leq C(\delta_{t}+h^{2}),$ (3.34)

and from (3.34) and (3.20) we have

$|\mathfrak{e}^{n}|_{1,h}\leq C(\delta_{t}+h^{2}).$ (3.35)

Then using (3.20), (3.28) and (3.35) in discrete Sobolev’s embedding theorem,
we derive that

$\parallel{e}^{n}\parallel_{\infty,h}+\parallel\mathfrak{e}^{n}\parallel_{\infty,h}\leq
C(\delta_{t}+h^{2}).$

Hence the theorem. ∎

######  Theorem 3.3 (Convergence of the scheme $\textbf{S}2$).

Let $y(x,t)$ and $p(x,t)$ be sufficiently smooth functions. For $\delta_{t}$
sufficiently small, the finite difference scheme (2.3)-(2.2) is first order in
time and second order in space convergent, i.e.,

$\max\limits_{n}\left\\{\parallel y^{n}-Y^{n}\parallel_{\infty,h}+\parallel
p^{n}-P^{n}\parallel_{\infty,h}\right\\}\leq C(\delta_{t}+h^{2}).$

###### Proof.

We provide an estimate that is consistent with the linear state approximation
scheme (2.3), and for an estimate that is consistent with discrete adjoint, we
turn to Theorem 3.2. Using Taylor expansion we can see that the exact solution
$y_{i}^{n}=y(x_{i},t_{n})$ and $p_{i}^{n}=p(x_{i},t_{n})$ satisfy the
following equations

$\frac{y_{i}^{n+1}-y_{i}^{n}}{\delta_{t}}=\Delta_{h}(y_{i}^{n})^{2}y_{i}^{n+1}-\Delta_{h}y_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}y_{i}^{n+1}+\frac{1}{\lambda}p_{i}^{n}+\widetilde{F}_{i}^{n},$
(3.36)

where $\widetilde{F}_{i}^{n}$ denotes the truncation error, which satisfies
the following for some positive constant $c_{3}$,

$\max\limits_{i,n}|\widetilde{F}_{i}^{n}|\leq c_{3}(\delta_{t}+h^{2}).$ (3.37)

Let us define the error $e_{i}^{n}=y_{i}^{n}-Y_{i}^{n}$ and
$\mathfrak{e}_{i}^{n}=p_{i}^{n}-P_{i}^{n}$. Then taking the difference between
(3.36) and (2.3) yields

$\frac{e_{i}^{n+1}-e_{i}^{n}}{\delta_{t}}=\Delta_{h}\left((y_{i}^{n})^{2}y_{i}^{n+1}-(Y_{i}^{n})^{2}Y_{i}^{n+1}\right)-\Delta_{h}e_{i}^{n}-\epsilon^{2}\Delta_{h}^{2}e_{i}^{n+1}+\frac{1}{\lambda}\mathfrak{e}_{i}^{n}+\widetilde{F}_{i}^{n}.$
(3.38)

Taking the inner product of (3.38) and $e^{n+1}$, we get

$\left(\frac{e^{n+1}-e^{n}}{\delta_{t}},e^{n+1}\right)_{h}=\left(\Delta_{h}\left((y^{n})^{2}y^{n+1}-(Y^{n})^{2}Y^{n+1}\right),e^{n+1}\right)-\left(\Delta_{h}e^{n},e^{n+1}\right)_{h}\\\
-\epsilon^{2}|e^{n+1}|_{2,h}^{2}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},e^{n+1}\right)_{h}+\left(\widetilde{F}^{n},e^{n+1}\right)_{h}.$
(3.39)

Observe that

$\left(\Delta_{h}\left((y^{n})^{2}y^{n+1}-(Y^{n})^{2}Y^{n+1}\right),e^{n+1}\right)\leq\left(|\widetilde{f}^{\prime}(y^{n})+\widetilde{f}^{\prime}(Y^{n})||e^{n+1}|,\Delta_{h}e^{n+1}\right)\leq
2\widetilde{M}\parallel e^{n+1}\parallel_{h}|e^{n+1}|_{2,h}.$

Using the above estimate, Young’s inequality and (3.37) on (3.39) we obtain

$\frac{\parallel e^{n+1}\parallel_{h}^{2}-\parallel
e^{n}\parallel_{h}^{2}}{2\delta_{t}}+\frac{\epsilon^{2}}{2}|e^{n+1}|_{2,h}^{2}\leq\underbrace{\left(\frac{4\widetilde{M}^{2}}{\epsilon^{2}}+\frac{1}{2\epsilon^{2}}+\frac{1}{2\lambda}+\frac{1}{2}\right)}_{\widetilde{\alpha}}\parallel
e^{n+1}\parallel_{h}^{2}\\\
+\frac{\epsilon^{2}}{2}|e^{n}|_{2,h}^{2}+\frac{1}{2\lambda}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{c_{3}}{2}(\delta_{t}+h^{2})^{2}.$
(3.40)

From (3.40) we have

$\left(1-2\delta_{t}\widetilde{\alpha}\right)\parallel
e^{n+1}\parallel_{h}^{2}+\epsilon^{2}\delta_{t}\left(|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}\right)\leq\parallel
e^{n}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+c_{3}\delta_{t}(\delta_{t}+h^{2})^{2}.$
(3.41)

For $\delta_{t}<\frac{1}{2\widetilde{\alpha}}$ in (3.41) and taking the sum
over $n$ we get

$\parallel
e^{n}\parallel_{h}^{2}\leq\frac{1}{\left(1-2\delta_{t}\widetilde{\alpha}\right)}\sum\limits_{i=1}^{n}\parallel
e^{i}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda\left(1-2\delta_{t}\widetilde{\alpha}\right)}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+\frac{c_{3}T}{\left(1-2\delta_{t}\widetilde{\alpha}\right)}(\delta_{t}+h^{2})^{2}.$
(3.42)

Applying Gronwall’s inequality on (3.42) for sufficiently small $\delta_{t}$
yields

$\parallel e^{n}\parallel_{h}^{2}\leq
C\delta_{t}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.43)

Adding (3.18) (the error estimates in Theorem 3.2 for adjoint equation) and
(3.43), and employing Gronwall’s inequality for sufficiently small
$\delta_{t}$ we have

$\parallel e^{n}\parallel_{h}+\parallel\mathfrak{e}^{n}\parallel_{h}\leq
C(\delta_{t}+h^{2}).$ (3.44)

Next we take the inner product of (3.38) and
$\partial_{t}e^{n}:=\frac{e^{n+1}-e^{n}}{\delta_{t}}$, that yields

$\displaystyle\parallel\partial_{t}e^{n}\parallel^{2}_{h}=$
$\displaystyle\left(\Delta_{h}\left((y^{n})^{2}y^{n+1}-(Y^{n})^{2}Y^{n+1}\right),\partial_{t}e^{n}\right)-\left(\Delta_{h}e^{n},\partial_{t}e^{n}\right)_{h}$
(3.45)
$\displaystyle-\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},\partial_{t}e^{n}\right)_{h}+\left(\widetilde{F}^{n},\partial_{t}e^{n}\right)_{h},$
$\displaystyle\leq 2\widetilde{M}\parallel
e^{n+1}\parallel_{h}\parallel\partial_{t}(\Delta_{h}e^{n})\parallel_{h}+|e^{n}|_{2,h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h}$
$\displaystyle-\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}+\frac{1}{2\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{1}{2}\parallel\partial_{t}e^{n}\parallel^{2}_{h}+\parallel\widetilde{F}^{n}\parallel_{h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h}.$

Using (3.37) and (3.20) on (3.45) we get

$\epsilon^{2}\frac{|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}}{2\delta_{t}}\leq\frac{\epsilon^{2}}{4}\frac{|e^{n+1}|_{2,h}^{2}+|e^{n}|_{2,h}^{2}}{\delta_{t}}+|e^{n}|_{2,h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.46)

From (3.46) we have

$\frac{\epsilon^{2}}{2}|e^{n+1}|_{2,h}^{2}\leq\left(\frac{3\epsilon^{2}}{2}+2\delta_{t}\right)|e^{n}|_{2,h}^{2}+2\delta_{t}C(\delta_{t}+h^{2})^{2}.$
(3.47)

Taking the sum over $n$ on (3.47) we get

$\frac{\epsilon^{2}}{2}|e^{n}|_{2,h}^{2}\leq\left(\frac{3\epsilon^{2}}{2}+2\delta_{t}\right)\sum\limits_{i=1}^{n}|e^{i}|_{2,h}^{2}+2TC(\delta_{t}+h^{2})^{2}.$
(3.48)

Applying Gronwall’s inequality on (3.48) we have

$|e^{n}|_{2,h}\leq C(\delta_{t}+h^{2}).$ (3.49)

Using (3.44) and (3.49) we obtain

$|e^{n}|_{1,h}\leq C(\delta_{t}+h^{2}).$ (3.50)

Now using (3.44), (3.50) and (3.35) in discrete Sobolev’s embedding theorem we
get our result. ∎

######  Theorem 3.4 (Convergence of the scheme $\textbf{S}3$).

Let $y(x,t)$ and $p(x,t)$ be sufficiently smooth functions. For $\delta_{t}$
sufficiently small, the finite difference scheme (2.4)-(2.2) is first order in
time and second order in space convergent, i.e.,

$\max\limits_{n}\left\\{\parallel y^{n}-Y^{n}\parallel_{\infty,h}+\parallel
p^{n}-P^{n}\parallel_{\infty,h}\right\\}\leq C(\delta_{t}+h^{2}).$

###### Proof.

We provide an estimate that is consistent with the linear state approximation
scheme (2.3), and for an estimate that is consistent with discrete adjoint, we
turn to Theorem 3.2. Using Taylor expansion we can see that the exact solution
$y_{i}^{n}=y(x_{i},t_{n})$ and $p_{i}^{n}=p(x_{i},t_{n})$ satisfy the
following equation

$\frac{y_{i}^{n+1}-y_{i}^{n}}{\delta_{t}}=\Delta_{h}(y_{i}^{n})^{3}-3\Delta_{h}y_{i}^{n}+2\Delta_{h}y_{i}^{n+1}-\epsilon^{2}\Delta_{h}^{2}y_{i}^{n+1}+\frac{1}{\lambda}p_{i}^{n}+\widehat{F}_{i}^{n},$
(3.51)

where $\widehat{F}_{i}^{n}$ denotes the truncation error, which satisfies the
following for some positive constant $c_{4}$

$\max\limits_{i,n}|\widehat{F}_{i}^{n}|\leq c_{4}(\delta_{t}+h^{2}).$ (3.52)

Let us define the error $e_{i}^{n}=y_{i}^{n}-Y_{i}^{n}$ and
$\mathfrak{e}_{i}^{n}=p_{i}^{n}-P_{i}^{n}$. Then taking the difference between
(3.51) and (2.4) yields

$\frac{e_{i}^{n+1}-e_{i}^{n}}{\delta_{t}}=\Delta_{h}\left((y_{i}^{n})^{3}-(Y_{i}^{n})^{3}\right)-3\Delta_{h}e_{i}^{n}+2\Delta_{h}e_{i}^{n+1}-\epsilon^{2}\Delta_{h}^{2}e_{i}^{n+1}+\frac{1}{\lambda}\mathfrak{e}_{i}^{n}+\widehat{F}_{i}^{n}.$
(3.53)

Taking the inner product of (3.53) and $e^{n+1}$ yields

$\displaystyle\left(\frac{e^{n+1}-e^{n}}{\delta_{t}},e^{n+1}\right)_{h}=$
$\displaystyle\left(\Delta_{h}\left((y_{i}^{n})^{3}-(Y_{i}^{n})^{3}\right),e^{n+1}\right)-3\left(\Delta_{h}e^{n},e^{n+1}\right)_{h}$
(3.54)
$\displaystyle+2\left(\Delta_{h}e^{n+1},e^{n+1}\right)_{h}-\epsilon^{2}|e^{n+1}|_{2,h}^{2}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},e^{n+1}\right)_{h}+\left(\widehat{F}^{n},e^{n+1}\right)_{h}.$

Using the differentiability of $f$ and Cauchy-Schwarz on (3.54) we have

$\displaystyle\frac{\parallel e^{n+1}\parallel_{h}^{2}-\parallel
e^{n}\parallel_{h}^{2}}{2\delta_{t}}$ $\displaystyle\leq M\parallel
e^{n}\parallel_{h}|e^{n+1}|_{2,h}+3\parallel
e^{n+1}\parallel_{h}|e^{n}|_{2,h}+2\parallel
e^{n+1}\parallel_{h}|e^{n+1}|_{2,h}$ (3.55)
$\displaystyle-{\epsilon^{2}}|e^{n+1}|_{2,h}^{2}+\frac{1}{\lambda}\parallel\mathfrak{e}^{n}\parallel_{h}\parallel
e^{n+1}\parallel_{h}+\parallel\widehat{F}^{n}\parallel_{h}\parallel
e^{n+1}\parallel_{h}$ $\displaystyle\leq\frac{M^{2}}{\epsilon^{2}}\parallel
e^{n}\parallel_{h}^{2}+\frac{\epsilon^{2}}{4}|e^{n+1}|_{2,h}+\frac{9}{2\epsilon^{2}}\parallel
e^{n+1}\parallel_{h}^{2}+\frac{\epsilon^{2}}{2}|e^{n}|_{2,h}+\frac{4}{\epsilon^{2}}\parallel
e^{n+1}\parallel_{h}^{2}+\frac{\epsilon^{2}}{4}|e^{n+1}|_{2,h}$
$\displaystyle-{\epsilon^{2}}|e^{n+1}|_{2,h}^{2}+\frac{1}{2\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{1}{2}\parallel
e^{n+1}\parallel_{h}^{2}+\frac{1}{2}\parallel\widehat{F}^{n}\parallel_{h}^{2}+\frac{1}{2}\parallel
e^{n+1}\parallel_{h}^{2},$
$\displaystyle\leq\frac{M^{2}}{\epsilon^{2}}\parallel
e^{n}\parallel_{h}^{2}+\underbrace{\left(\frac{17}{2\epsilon^{2}}+1\right)}_{\alpha_{0}}\parallel
e^{n+1}\parallel_{h}^{2}-\frac{\epsilon^{2}}{2}|e^{n+1}|_{2,h}^{2}+\frac{\epsilon^{2}}{2}|e^{n}|_{2,h}^{2}$
$\displaystyle\quad\quad\quad+\frac{1}{2\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{c_{4}}{2}(\delta_{t}+h^{2})^{2},$

where on the second inequality we use Young’s inequality and in the last
inequality we use (3.52). From (3.55) we get

$\left(1-2\delta_{t}\alpha_{0}\right)\parallel
e^{n+1}\parallel_{h}^{2}+\epsilon^{2}\delta_{t}\left(|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}\right)\leq\left(1+\frac{2\delta_{t}M^{2}}{\epsilon^{2}}\right)\parallel
e^{n}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+c_{4}\delta_{t}(\delta_{t}+h^{2})^{2}.$
(3.56)

For $\delta_{t}<\frac{1}{2{\alpha_{0}}}$ we take the sum over $n$ to obtain

$\parallel
e^{n}\parallel_{h}^{2}\leq\frac{\left(1+\frac{2\delta_{t}M^{2}}{\epsilon^{2}}\right)}{\left(1-2\delta_{t}{\alpha_{0}}\right)}\sum\limits_{i=1}^{n}\parallel
e^{i}\parallel_{h}^{2}+\frac{\delta_{t}}{\lambda^{2}\left(1-2\delta_{t}{\alpha_{0}}\right)}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+\frac{c_{4}T}{\left(1-2\delta_{t}{\alpha_{0}}\right)}(\delta_{t}+h^{2})^{2}.$
(3.57)

Applying Gronwall’s inequality on (3.57) for sufficiently small $\delta_{t}$,
we get

$\parallel e^{n}\parallel_{h}^{2}\leq
C\delta_{t}\sum\limits_{i=0}^{n-1}\parallel\mathfrak{e}^{i}\parallel_{h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.58)

Adding (3.18) (error estimates in Theorem 3.2 for adjoint equation) and (3.58)
, and employing Gronwall’s inequality for sufficiently small $\delta_{t}$ we
have

$\parallel e^{n}\parallel_{h}+\parallel\mathfrak{e}^{n}\parallel_{h}\leq
C(\delta_{t}+h^{2}).$ (3.59)

Next we take the inner product of (3.53) and
$\partial_{t}e^{n}:=\frac{e^{n+1}-e^{n}}{\delta_{t}}$, that produces

$\displaystyle\parallel\partial_{t}e^{n}\parallel^{2}_{h}=$
$\displaystyle\left(\Delta_{h}\left((y_{i}^{n})^{3}-(Y_{i}^{n})^{3}\right),\partial_{t}e^{n}\right)-3\left(\Delta_{h}e^{n},\partial_{t}e^{n}\right)_{h}+2\left(\Delta_{h}e^{n+1},\partial_{t}e^{n}\right)_{h}$
(3.60)
$\displaystyle-\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}+\frac{1}{\lambda}\left(\mathfrak{e}^{n},\partial_{t}e^{n}\right)_{h}+\left(\widehat{F}^{n},\partial_{t}e^{n}\right)_{h},$
$\displaystyle\leq{M}\parallel
e^{n}\parallel_{h}\parallel\partial_{t}(\Delta_{h}e^{n})\parallel_{h}+9|e^{n}|_{2,h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h}+4|e^{n+1}|_{2,h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h}$
$\displaystyle-\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}+\frac{1}{\lambda^{2}}\parallel\mathfrak{e}^{n}\parallel_{h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h}+\parallel\widehat{F}^{n}\parallel_{h}^{2}+\frac{1}{4}\parallel\partial_{t}e^{n}\parallel^{2}_{h},$

where on the second inequality we use differentiability of $f$ and Young’s
inequality. Using (3.59), (3.52) and (3.20) on (3.60) we have

$\epsilon^{2}\left(\Delta_{h}e^{n+1},\partial_{t}(\Delta_{h}e^{n})\right)_{h}\leq\frac{\epsilon^{2}}{4}\parallel\partial_{t}(\Delta_{h}e^{n})\parallel_{h}^{2}+9|e^{n}|_{2,h}^{2}+4|e^{n+1}|_{2,h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.61)

From (3.61) we get

$\epsilon^{2}\frac{|e^{n+1}|_{2,h}^{2}-|e^{n}|_{2,h}^{2}}{2\delta_{t}}\leq\frac{\epsilon^{2}}{4}\frac{|e^{n+1}|_{2,h}^{2}+|e^{n}|_{2,h}^{2}}{\delta_{t}}+9|e^{n}|_{2,h}^{2}+4|e^{n+1}|_{2,h}^{2}+C(\delta_{t}+h^{2})^{2}.$
(3.62)

This further implies

$\left(\frac{\epsilon^{2}}{2}-4\delta_{t}\right)|e^{n+1}|_{2,h}^{2}\leq\left(\frac{3\epsilon^{2}}{2}+9\delta_{t}\right)|e^{n}|_{2,h}^{2}+2\delta_{t}C(\delta_{t}+h^{2})^{2}.$
(3.63)

For $\delta_{t}<\frac{\epsilon^{2}}{8}$ we now take the sum over $n$ in (3.63)
to obtain

$\left(\frac{\epsilon^{2}}{2}-4\delta_{t}\right)|e^{n}|_{2,h}^{2}\leq\left(\frac{3\epsilon^{2}}{2}+9\delta_{t}\right)\sum\limits_{i=1}^{n}|e^{i}|_{2,h}^{2}+2TC(\delta_{t}+h^{2})^{2}.$
(3.64)

Applying Gronwall’s inequality on (3.64) we have

$|e^{n}|_{2,h}\leq C(\delta_{t}+h^{2}).$ (3.65)

Using (3.59) and (3.65) we obtain

$|e^{n}|_{1,h}\leq C(\delta_{t}+h^{2}).$ (3.66)

Then using (3.59), (3.66) and (3.35) in discrete Sobolev’s embedding theorem
we get our theorem. ∎

######  Remark 3.5.

To obtain the stability estimates of the proposed schemes in
$\parallel\cdot\parallel_{h}$, one can take an inner product of discrete state
equations (2.1), (2.3) and (2.4) with $Y^{n+1}$ and inner product of discrete
adjoint equation (2.2) with $P^{n}$. By following the procedures of the
aforementioned convergence proofs for all the schemes, one can obtain the
estimate for $\delta_{t}\leq O(\epsilon^{2})$, $\parallel
Y^{n}\parallel_{h}+\parallel P^{n}\parallel_{h}\leq
C\sum\limits_{i=1}^{n}\parallel\widehat{Y}^{i}\parallel_{h}$ for some generic
constant $C$. However, in numerical experiments we observe that even for
$\delta_{t}>O(\epsilon^{2})$ the schemes are stable and accurate.

######  Remark 3.6.

Derivation of the proposed schemes in higher dimensions is straightforward in
finite difference setting and proceeds in precisely the same manner as in the
1D case, so we omit the specifics. We only show the numerics in higher
dimensions.

######  Remark 3.7.

One can replace the term $f^{\prime}(Y_{i}^{n})\Delta_{h}P_{i}^{n}$ in (2.2)
by $f^{\prime}(Y_{i}^{n+1})\Delta_{h}P_{i}^{n}$ and get the same convergence
results for the respective schemes.

##  4\. Numerical Illustration

We now investigate the solution behaviour of the OCP (1.1)-(1.2) using the
proposed discretization strategies. To solve the nonlinear scheme
$\textbf{S}1$ we use Newton method at each time step with residual history of
$1e{-10}$.

###  4.1. Experiments in 1D

First we run experiments in 1D and demonstrate the accuracy and stability of
those proposed schemes in space and time. Moreover we run the schemes for
different desired state.

####  4.1.1. Accuracy test

Since the exact solution of the OCP (1.1)-(1.2) is not known, a comparison
between the solution of (1.1)-(1.2) on a coarse mesh with that on a fine mesh
is considered to check the accuracy. The error in the numerical solution is
measured in $L^{\infty}(0,T;L^{\infty}(\Omega))$. We take the spatial
computational domain $\Omega=(0,1)$. The initial and desired state for the
test are $y(x,0)=\cos(2\pi x)$ and $\widehat{y}(x,t)=\cos(2\pi x)e^{-t}$
respectively. First, we check the precision with respect to step size for all
the proposed schemes for a time window $[0,T=0.1]$. We discretize the spatial
domain by taking a mesh size $h=1/256$.

(a)

(b)

Figure 1. Convergence curves in time in $\log-\log$ scale with
$\epsilon=0.05,\lambda=0.1$ : on the left for the state equation, on the right
for the adjoint equation

In Figure 1 we observe that the first order convergence rate in time for both
state and adjoint equation is achieved for all the proposed schemes. Also note
that schemes are stable and accurate for $\delta_{t}>O(\epsilon^{2})$. To test
the accuracy in mesh size we fix the final time at $T=0.01$ and $N_{t}=200$.

(a)

(b)

Figure 2. Convergence curves in mesh size in $\log-\log$ scale with
$\epsilon=0.05,\lambda=0.1$ : on the left for the state equation, on the right
for the adjoint equation

Figure 2 shows how the suggested schemes accomplish the second order
convergence in mesh size for both state and adjoint equations.

####  4.1.2. Solution for different target states

In this subsection, we plot solutions of the state equation and control for
different desired state. For each experiment in this part, we use
$\Omega=(0,1)$ with mesh size $h=1/256$, time step $\delta_{t}=1e{-4}$ with
$N_{t}=100$. We run the scheme $\textbf{S}1$ to achieve the desired state
$\widehat{y}(x,t)=\cos(2\pi x)$.

(a)

(b)

(c)

Figure 3.  On the left : the desired state; In the middle: the computed state;
On the right : the computed control; For $\epsilon=0.05,\lambda=0.1$.

(a)

(b)

Figure 4.  On the left : the computed state; On the right : the computed
control; For $\epsilon=0.05,\lambda=1e{-4}$.

In Figure 3 we plot the desired state, the computed state and the optimal
control for $\epsilon=0.05$ and $\lambda=0.1$. To get an improvement in
accuracy on computed state we take $\lambda=1e{-4}$ and plot the solution in
Figure 4. We note that just the magnitude of calculated control changes, with
no improvement in computed state. Now, we test to check whether the calculated
state improves by changing $\epsilon$ to $\epsilon=0.09$. The solutions for
$\epsilon=0.09$ are shown in Figure 5, which demonstrates that the calculated
state is extremely close to the target state. So It is evident that in order
to get a solution profile that is identical to the intended state, the
calculated state relies on the choice of $\epsilon$.

(a)

(b)

Figure 5.  On the left : the computed state; On the right : the computed
control; For $\epsilon=0.09,\lambda=1e{-4}$.

Next we run the scheme $\textbf{S}2$ with the target state
$\widehat{y}(x,t)=\sin(\pi x)(t^{2}+1)$.

(a)

(b)

(c)

Figure 6.  On the left : the desired state; On the middle: the computed state;
On the right : the computed control; for $\epsilon=0.09,\lambda=1e{-4}$.

In Figure 6 we plot the desired state, the computed state and the optimal
control for $\epsilon=0.09$ and $\lambda=0.0001$. It is evident that the
calculated state is quite close to the target state. Lastly we run the final
scheme $\textbf{S}3$ with the desired state $\widehat{y}(x,t)=0.1e^{\cos(\pi
x)}(3t^{2}+1)$. Figure 7 shows the target state, calculated state, and optimal
control for $\epsilon=0.05$ and $\lambda=1e{-4}$, with an acceptable degree of
precision.

(a)

(b)

(c)

Figure 7.  On the left : the desired state; On the middle: the computed state;
On the right : the computed control; For $\epsilon=0.05,\lambda=1e{-4}$.

###  4.2. Experiments in 2D

For experiments in 2D we consider the spatial domain $\Omega=(0,1)^{2}$ with
uniform mesh size $h=1/50$ on both direction for all the experiments in this
subsection. First we run the scheme $\textbf{S}1$ with the desired state
$\widehat{y}(x,t)=\cos(2\pi xy)$ over the time window $[0,T=0.01]$ with step
size $\delta_{t}=5e{-4}$.

(a)

(b)

Figure 8.  On the left : the desired state; On the right : the computed state;
for $\epsilon=0.1,\lambda=1e{-4}$.

(a)

(b)

Figure 9.  On the left : the computed control at $t=0$; On the right : the
computed control at $t=T-\delta_{t}$; for $\epsilon=0.1,\lambda=1e{-4}$.

(a)

(b)

Figure 10.  On the left : the computed state at $t=T$; On the right : the
computed control at $t=T-\delta_{t}$; for $\epsilon=0.07,\lambda=1e{-4}$.

We illustrate the intended state and calculated state in Figure 8. We depict
the calculated control at $t=0$ and $t=T-\delta_{t}$ in Figure 9. We use
$\epsilon=0.07$ to enhance the calculated state, as seen in the left panel of
Figure 10; we additionally draw the corresponding control at $t=T-\delta_{t}$.

The scheme $\textbf{S}2$ is then executed to obtain the desired state
$\widehat{y}(x,t)=0.5\cos(2\pi x)\cos(2\pi y)(e^{-0.1t}+2t)$ over the time
window $[0,T=0.03]$ with time step $\delta_{t}=1e{-3}$

(a)

(b)

Figure 11.  On the left : the desired state; On the right : the computed
state; for $\epsilon=0.1,\lambda=0.01$.

(a)

(b)

Figure 12.  On the left : computed control at $t=0$; On the right : computed
control at $t=T-\delta_{t}$; for $\epsilon=0.1,\lambda=0.01$. (a) Figure 13.
Difference between state and control at $t=T$ for $\epsilon=0.1,\lambda=0.01$.

The final target state and calculated state are shown in Figure 11, and it is
evident from the plot that the two plots can not be distinguished. At $t=0$
and $t=T-\delta_{t}$ in Figure 212 we display the corresponding control. In
Figure 13, we exhibit the difference between the two states at the final time
to show how near the calculated state is to the desired state. The difference,
as we can see, is in the order of $10^{-2}$. Keep in mind that by selecting
the relatively larger values of $\epsilon=0.1$ and $\lambda=0.01$, we are able
to achieve the target state with an absolute difference of order $10^{-2}$. At
this point, we may thus conclude that, in addition to the choices of
$\epsilon$ and $\lambda$, the target state choice also has an impact.

The scheme $\textbf{S}3$ is now executed to attain the desired state
$\widehat{y}(x,t)=0.5\sin(xy)(1-2t)$ across the time frame $[0,T=0.01]$ with a
step sige of $\delta_{t}=3.33e{-4}$.

(a)

(b)

Figure 14.  On the left : the desired state; On the right : the computed
state; for $\epsilon=0.08,\lambda=0.001$.

(a)

(b)

Figure 15.  On the left : computed control at $t=0$; On the right : computed
control at $t=T-\delta_{t}$; for $\epsilon=0.08,\lambda=0.001$.

The final target state and calculated state at final time are shown in Figure
14. It is obvious from the plot that the computed state is close to the
desired state. We show the corresponding control at $t=0$ and $t=T-\delta_{t}$
in Figure 15.

##  5\. Conclusions

In this work, we formulate linear and nonlinear discretization schemes for the
OCP with the CH equation as constrain. We present rigorous convergence
analysis for all the proposed schemes. Lastly we verify the numerical accuracy
and present numerical solution for various desired state in 1D and 2D.

## Acknowledgement

The authors would like to thank the CSIR (File No : 09/1059(0019)/2018-EMR-I)
and DST-SERB (File No : SRG/2019/002164) for the research grant and IIT
Bhubaneswar for providing excellent research environment.

## References

  * [1] J. W. Cahn, On spinodal decomposition, Acta Metall, 9 (1961), pp. 795–801.
  * [2] J. W. Cahn and W. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, J. Chem. Phys., 28 (1958), pp. 258–267.
  * [3] K. Cheng, W. Feng, C. Wang, and S. M. Wise, An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation, J. Comput. Appl. Math., 362 (2019), pp. 574–595.
  * [4] A. Christlieb, K. Promislow, and Z. Xu, On the unconditionally gradient stable scheme for the Cahn-Hilliard equation and its implementation with Fourier method, Commun. Math. Sci., 11 (2013), pp. 345–360.
  * [5] Q. Du and R. A. Nicolaides, Numerical analysis of a continuum model of phase transition, SIAM J. Numer. Anal., 28 (1991), pp. 1310–1322.
  * [6] N. Duan and X. Zhao, Optimal control for the multi-dimensional viscous cahn–hilliard equation, Electron. J. Differential Equations, (2015), p. 13.
  * [7] C. M. Elliott, The Cahn-Hilliard model for the kinetics of phase separation, in Mathematical models for phase change problems (Óbidos, 1988), vol. 88 of Internat. Ser. Numer. Math., Birkhäuser, Basel, 1989, pp. 35–73.
  * [8] C. M. Elliott and D. A. French, Numerical studies of the cahn-hilliard equation for phase separation, IMA Journal of Applied Mathematics, 38 (1987), pp. 97–128.
  * [9] C. M. Elliott and Z. Songmu, On the Cahn-Hilliard equation, Arch. Rational Mech. Anal., 96 (1986), pp. 339–357.
  * [10] D. J. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation, in Computational and mathematical models of microstructural evolution (San Francisco, CA, 1998), vol. 529 of Mater. Res. Soc. Sympos. Proc., MRS, Warrendale, PA, 1998, pp. 39–46.
  * [11] D. Furihata, A stable and conservative finite difference scheme for the cahn-hilliard equation, Numerische Mathematik, 87 (2001), pp. 675–699.
  * [12] J. G. Heywood and R. Rannacher, Finite-element approximation of the nonstationary navier–stokes problem. part iv: error analysis for second-order time discretization, SIAM Journal on Numerical Analysis, 27 (1990), pp. 353–384.
  * [13] M. Hintermüller and D. Wegner, Distributed optimal control of the cahn–hilliard system including the case of a double-obstacle homogeneous free energy density, SIAM Journal on Control and Optimization, 50 (2012), pp. 388–418.
  * [14] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, Optimization with PDE constraints, vol. 23, Springer Science & Business Media, 2008.
  * [15] D. Lee, J.-Y. Huh, D. Jeong, J. Shin, A. Yun, and J. Kim, Physical, mathematical, and numerical derivations of the cahn–hilliard equation, Computational Materials Science, 81 (2014), pp. 216–225.
  * [16] J. Shin, D. Jeong, and J. Kim, A conservative numerical method for the cahn–hilliard equation in complex domains, Journal of Computational Physics, 230 (2011), pp. 7441–7455.
  * [17] A. M. Stuart and A. R. Humphries, Model problems in numerical stability theory for initial value problems, SIAM Rev., 36 (1994), pp. 226–257.
  * [18] F. Tröltzsch, Optimal control of partial differential equations: theory, methods, and applications, vol. 112, American Mathematical Soc., 2010\.
  * [19] J. Yong and S. Zheng, Feedback stabilization and optimal control for the cahn-hilliard equation, Nonlinear Analysis: Theory, Methods & Applications, 17 (1991), pp. 431–444.
  * [20] X. Zhao and C. Liu, Optimal control problem for viscous cahn–hilliard equation, Nonlinear Analysis: Theory, Methods & Applications, 74 (2011), pp. 6348–6357.
  * [21] , Optimal control for the convective cahn–hilliard equation in 2d case, Applied Mathematics & Optimization, 70 (2014), pp. 61–82.
  * [22] J. Zheng and Y. Wang, Optimal control problem for cahn–hilliard equations with state constraint, Journal of Dynamical and Control Systems, 21 (2015), pp. 257–272.