Optimal control of multiphase free boundary problems for nonlinear parabolic equations

TL;DR

This paper develops a method for optimal boundary heat flux control in multiphase free boundary problems modeled by nonlinear parabolic PDEs, ensuring convergence of discrete approximations to the continuous problem.

Contribution

It introduces a finite difference approach for controlling multiphase Stefan problems and proves the existence and convergence of optimal controls.

Findings

Established existence of optimal control.

Proved convergence of discrete controls to continuous control.

Derived uniform bounds and energy estimates for the PDE solutions.

Abstract

We consider the optimal control of singular nonlinear partial differential equation which is the distributional formulation of the multiphase Stefan type free boundary problem for the general second order parabolic equation. Boundary heat flux is the control parameter, and the optimality criteria consist of the minimization of the $L_2$-norm declination of the trace of the solution to the PDE problem at the final moment from the given measurement. Sequence of finite-dimensional optimal control problems is introduced through finite differences. We establish existence of the optimal control and prove the convergence of the sequence of discrete optimal control problems to the original problem both with respect to functional and control. Proofs rely on establishing a uniform $L_{\infty}$ bound, and $W_2^{1,1}$-energy estimate for the discrete nonlinear PDE problem with discontinuous coefficient.

Equations

$$\mathcal{J}(g)=\Vert v(x,T;g)-\omega (x){\Vert }_{L_2(0,\ell )}^2$$

$$\mathscr{G}_{R}=\big{\{}g:g\in W_{2}^{1}(0,T),\,\,\|g\|_{W_{2}^{1}[0,T]}\leq R\big{\}}.$$

$$\frac{\partial \beta (v)}{\partial t}-\mathcal{L}v-f(x,t)\ni 0,$$

$$v(x,0)=\Phi (x),$$

$$a{v}_x+bv{|}_{x=0}=g(t),$$

$$a{v}_x+bv{|}_{x=\ell }=p(t),$$

$$\beta(r)=\left\{\begin{array}[]{l}\beta_{j}(r)+\sum\limits_{i=0}^{j-1}\nu_{i},\quad\text{for}\ v^{j-1}<r<v^{j},\\ $\Big{[}$\beta_{j}(v^{j})+\sum\limits_{i=0}^{j-1}\nu_{i},\beta_{j}(v^{j})+\sum\limits_{i=0}^{j}\nu_{i}$\Big{]}$,\quad\text{for}\ \ r=v^{j},\\ \beta_{j+1}(r)+\sum\limits_{i=0}^{j}\nu_{i},\quad\text{for}\ v^{j}<r<v^{j+1};\ j=1,2,...,m\end{array}\right.$$

$$0<\bar{b}\le {\beta }_j^{\prime }(r)\le \bar{B},j=1,...,m+1;$$

$$\mathcal{L}v=(a(x,t)v_x+b(x,t)v{)}_x-c(x,t)v$$

$$a(x,t)\ge a_0>0,\text{a.e.}(x,t)\in D=\{0<x<\ell ,0<t<T\}.$$

$$(u,v)={\int }_0^{T}uvdt.$$

$$(u,v)={\int }_{D}uvdxdt.$$

$$\Vert u{\Vert }_{L_{\infty }[0,T]}=\underset{0\le t\le T}{\mathrm{esssup}}|u(t)|.$$

$$\Vert u{\Vert }_{L_{\infty }(D)}=\underset{(x,t)\in D}{\mathrm{esssup}}|u(x,t)|.$$

$$(u,v)={\int }_0^T\sum _{s=0}^k\frac{d^{s}u}{d{t}^s}\frac{d^{s}v}{d{t}^s}dt.$$

$$(u,v)=\int_{D}\Big{(}uv+\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}\Big{)}\,dx\,dt.$$

$$(u,v)=\int_{D}\Big{(}uv+\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial t}\frac{\partial v}{\partial t}\Big{)}\,dx\,dt.$$

$$\Vert u{\Vert }_{W_{\infty }^{1,0}(D)}=\Vert u{\Vert }_{L_{\infty }(D)}+\Vert u_x{\Vert }_{L_{\infty }(D)}$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{0}^{\ell}\Big{[}-B(x,t,v(x,t))\psi_{t}+a(x,t)v_{x}\psi_{x}+b(x,t)v\psi_{x}+c(x,t)v\psi-f\psi\Big{]}\,dxdt$$

$$-\underset{0}{\overset{\ell }{\int }}{B}_0(x,0,\Phi (x))\psi (x,0)dx-\underset{0}{\overset{T}{\int }}p(t)\psi (\ell ,t)dt+\underset{0}{\overset{T}{\int }}g(t)\psi (0,t)dt=0,$$

$${\omega }_{\tau }=\{t_k,k=\overline{1,n}\},\tau =\frac{T}{n},t_k=k\tau ,{\omega }_h=\{x_i,i=\overline{1,m}\},h=\frac{\ell }{m},x_i=ih$$

$$\frac{h}{\tau }\ge \frac{8\Vert b{\Vert }_{L_{\infty }}}{\bar{b}}$$

$$w_k=\frac{1}{\tau }\underset{t_{k-1}}{\overset{t_k}{\int }}w(t)dt,{\Phi }_i=\frac{1}{h}\underset{x_i}{\overset{x_{i+1}}{\int }}\Phi (x)dx,\text{leavevmode}\text{leavevmode}{\Phi }_m=\Phi (\ell ),$$

$$q_{ik}=\frac{1}{\tau h}\underset{t_{k-1}}{\overset{t_k}{\int }}\underset{x_i}{\overset{x_{i+1}}{\int }}q(x,t)dxdt,k=\overline{1,n},i=\overline{0,m-1},$$

$${\mathcal{G}}_R^n=\{[g{]}_n\in {\mathbb{R}}^{n+1}:\Vert [g{]}_n{\Vert }_{w_2^1}\le R\}$$

$$\Vert [g{]}_n{\Vert }_{w_2^1}^2=\sum _{k=1}^n\tau g_k^2+\sum _{k=1}^n\tau g_{k\bar{t}}^2$$

$${\mathcal{Q}}_n(g)=[g{]}_n,\text{for }g\in {\mathcal{G}}_R,\text{where }{g}_k=\frac{1}{\tau }\underset{t_{k-1}}{\overset{t_k}{\int }}g(t)dt,k=\overline{0,n},$$

$${\mathcal{P}}_n([g{]}_n)=g^n,\text{leavevmode}\text{for }[g{]}_n\in {\mathcal{G}}_R^n;$$

$$g^n(t)=g_{k-1}+\frac{g_k-g_{k-1}}{\tau }(t-t_{k-1}),t\in [t_{k-1},t_k),k=\overline{1,n}.$$

$$b_n(v)={\int }_{v-\frac{1}{n}}^{v+\frac{1}{n}}\beta (y){\omega }_n(v-y)dy,$$

$${\omega }_n(v)=\{\begin{array}{cc}\mathcal{C}n{e}^{-\frac{1}{1-n^2{v}^2}},& |v|\le \frac{1}{n}\\ 0,& |v|>\frac{1}{n}\end{array}$$

$$b_n^{\prime }(v)={\int }_{v-\frac{1}{n}}^{v+\frac{1}{n}}{\beta }^{\prime }(y){\omega }_n(v-y)dy.$$

$$b_n^{\prime }(v)\ge \bar{b}>0$$

$$\displaystyle\sum\limits_{i=0}^{m-1}h\Big{[}\big{(}b_{n}(v_{i}(k))\big{)}_{\bar{t}}\eta_{i}+a_{ik}v_{ix}(k)\eta_{ix}+b_{ik}v_{i}(k)\eta_{ix}+c_{ik}v_{i}(k)\eta_{i}-f_{ik}\eta_{i}\Big{]}$$

Full text

1. Introduction

1.1. Optimal Control Problem

Let $v^{1}<v^{2}<\cdots<v^{m}$ are given real numbers. Consider an optimal control problem on the minimization of the cost functional

[TABLE]

over the control set:

[TABLE]

where $v$ is a solution of the singular nonlinear PDE problem

[TABLE]

with $\beta(\cdot)$ being a maximal monotone graph of the form

[TABLE]

with a given positive constants $\nu_{j},j=1,...,m$; $\nu_{0}=0$, $v^{0}=-\infty$, $v^{m+1}=+\infty$; $\beta_{i}(\cdot),i=1,...,m+1$ are monotone increasing Lipschitzian functions in their respective domain of definition, $\beta_{j}(v^{j})=\beta_{j+1}(v^{j}),j=1,...,m$,

[TABLE]

and $\mathcal{L}$ is an elliptic operator

[TABLE]

with bounded measurable coefficients $a,b,c$ and

[TABLE]

Described optimal control problem will be called a Problem $\mathcal{S}$.

In the particular case of $\mathcal{L}=\Delta$, system (1.2)-(1.5) presents distributional formulation of the multiphase Stefan problem describing flow of the heat in the presence of phase transitions [41, 33, 35]. In the physical context, $v(x,t)$ is a temperature, $f(x,t)$ is a density of heat sources, $\Phi(x)$ is an initial temperature distribution, $g(t)$ and $p(t)$ are heat flux on fixed boundaries, $v^{j}$’s are phase transition temperatures; $\beta_{j}^{\prime}(v),v^{j}<v<v^{j+1},j=0,1,...,m$ characterize heat conductivities in each phase, and the positive jump constants $\nu^{j},j=1,...,m$ are expressing latent heat of fusion during phase transition. In particular, the case $m=1,v^{1}=0$ is a classical two-phase Stefan problem describing melting of the ice or freezing of the water. [25, 37]. More complex examples of multiphase Stefan problem includes biomedical problem about the laser ablation of biomedical tissues, which motivates general elliptic operator $\mathcal{L}$ with coefficients $a,b,c$ expressing anisotropic properties of the media. Optimal control problem (1.1)-(1.5) aims to achieve the desired temperature distribution $\omega(x)$ at the final moment by controlling the boundary flux $g(t)$ on the fixed boundary. Equivalently, it is a variational formulation of the inverse multiphase Stefan type free boundary problem on the identification boundary flux $g(t)$ through measurement $\omega(x)$ of the final moment temperature distribution.

The goal of this paper is to prove the well-posedness of the optimal control problem (1.1)-(1.5), and to prove the convergence of the sequence of the finite-dimensional discretized optimal control problems to the optimal control problem (1.1)-(1.5) both with respect to functional and control via the method of finite differences.

The idea of transformation of the multiphase Stefan problem to boundary value problem for singular PDE (1.2) originated in [41]. Existence and uniqueness of the weak solution was proved in [41, 33, 35] when $\mathcal{L}=\Delta$. Hölder continuity of the weak solutions was proved in [21, 20] for general nonlinear elliptic operators $\mathcal{L}$. Continuity of the weak solution for the two-phase Stefan problem was proved in [16].

The one-phase inverse Stefan problem (ISP) was first mentioned in [18], where phase transition boundary is known and heat flux on left boundary is to be found, and the variational approach for solving the ISP was used in [14, 15]. In [46] ISP was formulated as an optimal control problem and the existence of the optimal control is proved. In [48], the Frechet derivative was found, the convergence of finite difference schemes was proved, and Tikhonov regularization was suggested in order to improve results. The following works on the ISP split into two different directions: ISPs with given fixed phase transition boundaries ([9, 11, 13, 17, 19, 22, 23, 28, 44, 26]), and ISPs with unknown phase transition boundaries ([10, 24, 29, 30, 31, 32, 34, 36, 40, 38, 42, 43, 45, 26]). We can refer to [26] for a full list of references for both types stated above, which include both linear and quasilinear parabolic equations.

In [1, 2] a new variational formulation of the the one-phase ISP was developed, in which optimal control framework was implemented, where the phase transition boundary is included in the control set along with the boundary heat flux. The sum of the $L_{2}$-norm declinations are minimized against the available measurements of temperature on the fixed boundary, available measurements of the free boundary location, and temperature at final moment. Important advantage of the new control theoretic approach is that it can handle situations where the phase transition temperature is not known explicitly, and is available only through measurement with possible error. Another major advantage of the new variational method suggested in [1, 2] is based on the fact that for a given control vector corresponding state vector solves PDE problem in a fixed region instead of full free boundary problem. This allows to reduce significantly computational cost of iterative numerical methods based on gradient type methods in Sobolev spaces. In [3, 4], Frechet differentiability in Sobolev-Besov spaces was proved and the formula for the Frechet gradient and optimality condition are derived. In [6, 8] gradient method was implemented in Hilbert-Besov spaces framework for the numerical solution of the ISP.

The new method developed in [1, 2] is not applicable to inverse multiphase free boundary problems. The reason is that the Stefan condition on the free boundary includes the saltus of the boundary flux from neighbouring phases, and by fixing free boundary as a control parameter the Stefan condition does not become a Neumann or Robin type boundary condition for the PDE. In a recent paper [5], a new variational method was introduced for the solution of the inverse multiphase Stefan problem. The IMSP is reformulated in a new optimal control framework in which the boundary is fixed, yet the state vector satisfies a nonlinear PDE with coefficients possessing jump discontinuities along phase transition boundaries. In [5] existence of the optimal control and convergence of the sequence of discretized optimal control problems via method of finite differences is proved. In [7], this framework is extended to the multidimensional IMSP.

1.2. Notation of Function Spaces

$L_{2}(0,T)$ - Space of Lebesgue square-integrable functions. It is a Hilbert space with inner product

[TABLE]

$L_{2}(D)$ - Hilbert space with inner product

[TABLE]

$L_{\infty}(0,T)$ - Space of essentially bounded functions. It is a Banach space with norm

[TABLE]

$L_{\infty}(D)$ - Space of essentially bounded functions. It is a Banach space with norm

[TABLE]

$W_{2}^{k}(0,T),k=1,2,...$ - Hilbert space of all elements of $L_{2}(0,T)$ whose weak derivatives up to order $k$ exist and belong to $L_{2}(0,T)$. The inner product is defined as

[TABLE]

$W_{2}^{1,0}(D)$ - Hilbert space of all elements of $L_{2}(D)$ that have a weak derivative in the $x$ direction, $\frac{\partial u}{\partial x}$, and such that it belongs to $L_{2}(D)$. The inner product is defined as

[TABLE]

$W_{2}^{1,1}(D)$ - Hilbert space of all elements of $L_{2}(D)$ with weak derivatives of first order, $\frac{\partial u}{\partial x}$, $\frac{\partial u}{\partial t}$. Also its weak derivatives must belong to $L_{2}(D)$. The inner product is defined as

[TABLE]

$W_{\infty}^{1,0}(D)$ - Space of elements of $L_{\infty}(D)$ with weak derivative in the $x$ direction existing and belonging to $L_{\infty}(D)$. It is a Banach space with norm

[TABLE]

1.3. Weak Solution of the Multiphase Free Boundary Problem

We now formulate the notion of the weak solution of the nonlinear multiphase parabolic free boundary problem (1.2)-(1.5).

Definition 1.1**.**

We say that a measurable function $B(x,t,v)$ is of type $\mathscr{B}$ if

1. (a)

$B(x,t,v)=\beta(v),\qquad v\neq v^{j},\quad\forall j=\overline{1,J}$ 2. (b)

$B(x,t,v)\in[\beta(v^{j})^{-},\beta(v^{j})^{+}],\qquad v=v^{j}$ for some $j$.

Definition 1.2**.**

$v\in W_{2}^{1,1}(D)\cap L_{\infty}(D)$ is called a weak solution of the problem (1.2)-(1.5) if for any two functions $B,B_{0}$ of type $\mathscr{B}$, the following integral identity is satisfied

[TABLE]

for all $\psi\in W_{2}^{1,1}(D)$ with $\psi(x,T)=0$.

1.4. Discrete Optimal Control Problem

Let

[TABLE]

be grids in the time and space domains, respectively, under the assumptions that $m\rightarrow\infty\quad$ as $n\rightarrow\infty$ and

[TABLE]

Define the Steklov averages

[TABLE]

where $w$ represents any of the functions $p$, $\Gamma$, $g$, or $g^{n}$, and $q$ represents any of the functions $a,b,c,$ and $f$. Introduce the discretized control set

[TABLE]

where $[g]_{n}=(g_{0},g_{1},\ldots,g_{n})$, and

[TABLE]

with $g_{k\bar{t}}=\frac{g_{k}-g_{k-1}}{\tau}.$ Assume that every element $g\in W_{2}^{1}(0,T)$ is extended on the interval $[-\tau,0]$ as a constant $g(0)$. Consider now the mappings between the discrete and continuous control sets, $\mathscr{Q}_{n}:W_{2}^{1}(0,T)\rightarrow\mathbb{R}^{n+1},\quad\mathscr{P}_{n}:\mathbb{R}^{n+1}\rightarrow W_{2}^{1}(0,T)$ as

[TABLE]

Approximate the function $\beta(v)$ by the infinitely differentiable sequence

[TABLE]

where $\omega_{n}$ is a standard mollifier defined as

[TABLE]

and the constant $\mathcal{C}$ is chosen so that $\int\limits_{\mathbb{R}}\omega_{1}(u)\,du=1$. Since $\beta^{\prime}(v)$ is piecewise-continuous, we also have

[TABLE]

This implies $b_{n}$ is also strict monotonically increasing and by (1.6) we have

[TABLE]

We now define a solution to the problem (1.2)-(1.5) in the discrete sense

Discrete State Vector. Given $[g]_{n}$, the vector function [v([g]_{n})]_{n}\\ =\big{(}v(0),v(1),\ldots,v(n)\big{)}; $\quad v(k)\in\mathbb{R}^{m+1},\quad k=0,\ldots,n$ is called a discrete state vector if

1. (a)

$v_{i}(0)=\Phi_{i},\quad i=\overline{0,m},$ 2. (b)

For arbitrary $k=1,\ldots,n$, the vector $v(k)\in\mathbb{R}^{m+1}$ satisfies

[TABLE]

Given $[g]_{n}\in\mathscr{G}^{n}_{R}$, the discrete cost functional $\mathscr{I}_{n}$ is defined as

[TABLE]

where $v_{i}(k)$ are components of the discrete state vector $[v([g]_{n})]_{n}$. Finite-dimensional optimal control problem on the minimization of $\mathscr{I}_{n}([g]_{n})$ on a control set $\mathscr{G}_{R}^{n}$ will be called Problem $\mathcal{S}_{n}$. We define

[TABLE]

Furthermore, the following interpolations will be considered:

[TABLE]

2. Main Results

Unless stated otherwise, throughout the paper we assume the following conditions are satisfied by the data:

[TABLE]

the coefficient $a$ satisfies (1.7); and $\Phi(x)=v^{j},j=1,...,m$ on a set of measure 0 in the interval $(0,\ell)$.

Theorem 2.1**.**

The optimal control problem $\mathcal{S}$ has a solution, that is, the set

[TABLE]

is not empty.

Theorem 2.2**.**

The sequence of discrete optimal control problems $\mathcal{S}_{n}$ approximates the optimal control problem $\mathcal{S}$ with respect to functional, that is,

[TABLE]

where

[TABLE]

If $[g]_{n_{\varepsilon}}\in\mathscr{G}_{R}^{n}$ is chosen such that

[TABLE]

then the sequence $g^{n}=\mathscr{P}_{n}([g]_{n_{\varepsilon}})$ has a subsequence convergent to some element $g_{*}\in\mathscr{G}_{*}$ weakly in $W_{2}^{1}(0,T)$ and strongly in $L_{2}(0,T)$. Moreover, the piecewise linear interpolations $\hat{v}^{\tau}$ of the corresponding discrete state vectors $[v([g]_{n_{\varepsilon}})]_{n}$ converge to the weak solution $v(x,t;g_{*})\in W_{2}^{1,1}(D)\cap L_{\infty}(D)$ of the singular PDE problem (1.2)-(1.5) weakly in $W_{2}^{1,1}(D)$, strongly in $L_{2}(D)$, and almost everywhere on $D$.

3. Preliminary Results

Lemma 3.1**.**

Given any $[g]_{n}\in\mathscr{G}^{n}$, and any $h,\tau$, a discrete state vector exists uniquely.

Proof. First we prove uniqueness by induction. For a given $[g]_{n}$, suppose $v$ and $\tilde{v}$ both are discrete state vectors. Due to definition of how a discrete state vector is constructed, we have that $v(0)=\tilde{v}(0)$. Now suppose that $v(k-1)=\tilde{v}(k-1)$ for some fixed $k\geq 1$. Since $v$ and $\tilde{v}$ both satisfy (1.17), subtract the identities for both $v$ and $\tilde{v}$, choosing $\eta=v(k)-\tilde{v}(k)$ to get:

[TABLE]

From here we get using Cauchy inequality with $\epsilon=a_{0}$:

[TABLE]

Absorbing to left hand side, and by using (1.16), we get:

[TABLE]

The whole summand is non-negative for sufficiently small $\tau$. Therefore, it is equal to 0, which implies that $v_{i}(k)=\tilde{v}_{i}(k),\quad\forall i=\overline{0,m}$. Hence, by induction, $v=\tilde{v}$.

Now we seek to prove existence through induction. Construct $v(0)$ through definition of a Discrete State Vector. Note that $v(0)$ is bounded since $\|v(0)\|\leq\|\Phi\|_{L_{\infty}[0,\ell]}$. Fix $k\geq 1$, and assume that $v(k-1)$ has been constructed so that (1.17) is satisfied for all $K<k$. Moreover, assume that each element of $v(k-1)$ is bounded. Through manipulation, the summation identity (1.17) is equivalent to solving the following system of non-linear equations:

[TABLE]

We will construct $v(k)$ by the method of successive approximations. Fix $h$ and $\tau$, and choose $v^{0}=v(k-1)$. Having obtained $v^{N}$, we search $v^{N+1}$ as a solution of the following:

[TABLE]

We now proceed to prove that the sequence $\{v^{N}\}$ converges to the unique solution of (3.4). Subtract (3.5) for $N$ and $N-1$ to get

[TABLE]

which can be transformed to

[TABLE]

where

[TABLE]

Due to (1.16), we have $\zeta_{n,N}^{i}\geq\bar{b},\leavevmode\nobreak\ \leavevmode\nobreak\ i=\overline{0,m-1}$. Let

[TABLE]

From (3.8), taking the first equation into consideration, we have:

[TABLE]

We have

[TABLE]

and

[TABLE]

by (1.9) and for sufficiently small $h$ and $\tau$. Thus, $0<\delta_{0}^{-1}<1$. Similarly,

[TABLE]

Through similar argument as with $\delta_{0}$, we derive that $0<\delta_{i}^{-1}<1$ for $i=\overline{1,m-1}$ for $0<h<<1$. For $i=m$, we get

[TABLE]

For $0<h<<1$, we can see that the term in left brackets will be close to 1, and due to (1.9), as before we derive that $0<\delta_{m}^{-1}<1$ for sufficiently small $h$ and $\tau$. Let $\delta=\max\limits_{i=\overline{1,m}}\delta_{i}^{-1}$. We thus have $\delta<1$, and

[TABLE]

Following the proof given in [5] (Lemma 1, Section 2) it follows that there exist finite limits

[TABLE]

Passing to limit as $N\to+\infty$ in (3.5), we derive that $v_{i}(k),i=\overline{1,m}$ is a unique solution of (3.4). $\square$

Given the existence and uniqueness of the discrete state vector for fixed $n$, we can uniquely define for each $k=1,\ldots,n$ the vector $\zeta_{k}$ whose $m$ components $\zeta_{k}^{i}$ are given by

[TABLE]

The following is a well known necessary and sufficient condition for the convergence of the discrete optimal control problems to continuous optimal control problem.

Lemma 3.2**.**

[46]** The sequence of discrete optimal control problems approximates the continuous optimal control problem if and only if the following conditions are satisfied:

• •

for arbitrary sufficiently small $\varepsilon>0$ there exists $M_{1}=M_{1}(\varepsilon)$ such that $\mathscr{Q}_{M}(g)\in\mathscr{G}^{M}_{R}$ for all $g\in\mathscr{G}_{R-\varepsilon}$ and $M\geq M_{1}$; and for any fixed $\varepsilon>0$ and for all $g\in\mathscr{G}_{R-\varepsilon}$ the following inequality is satisfied:

[TABLE]

• •

for arbitrary sufficiently small $\varepsilon>0$ there exists $M_{2}=M_{2}(\varepsilon)$ such that $\mathscr{P}_{M}([g]_{M})\in\mathscr{G}_{R+\varepsilon}$ for all $[g]_{M}\in\mathscr{G}^{M}_{R}$ and $M\geq M_{2}$; and for all $[g]_{M}\in\mathscr{G}^{M}_{R}$, $M\geq 1$ the following inequality is satisfied:

[TABLE]

• •

the following inequalities are satisfied:

[TABLE]

where $\mathscr{J}_{*}(\pm\varepsilon)=\inf\limits_{\mathscr{G}_{R\pm\varepsilon}}\mathscr{J}(g)$.

Lemma 3.3**.**

[5]** The mappings $\mathscr{P}_{n},\mathscr{Q}_{n}$ satisfy the conditions of Lemma 3.2.

Lemma 3.4**.**

There is at most one solution to the multiphase free boundary problem (1.2)-(1.5) in the sense of (1.8).

Proof. The uniqueness of the weak solution is proved in Section $9$ of Chapter V of [35] for the classical multiphase Stefan Problem ($\mathcal{L}=\Delta$) under zero Dirichlet boundary conditions on the fixed boundary. Lemma 4 of [5] generalized the result to the case of non-homogeneous Neumann boundary condition. We generalize the result to the case of multiphase free boundary problem with general elliptic operator $\mathcal{L}$. Uniqueness is proved over a wider class of solutions than given in (1.8). Suppose that $v\in L_{\infty}(D)$ only, not necessarily in the Sobolev space $W_{2}^{1,1}(D)$, and that for any two functions $B,B_{0}$ of type $\mathscr{B}$ it satisfies the identity

[TABLE]

Any function satisfying (1.2) will also satisfy the above definition. Suppose $v$ and $\tilde{v}$ are two solutions in the sense of (3.17), and subtract (3.17) with solution $\tilde{v}$ from that of $v$. Due to $\Phi$ taking on phase transition temperatures on sets of measure 0, the $B_{0}$ term will vanish and we are left with the following:

[TABLE]

where $z(x,t)=\frac{v-\tilde{v}}{B(x,t,v)-\tilde{B}(x,t,\tilde{v})}$. For $(x,t)\in D$ such that $v(x,t)=\tilde{v}(x,t)$, we have $z(x,t)=0$. Otherwise, since $B$ and $\tilde{B}$ are strictly increasing on $v$ a.e. $(x,t)\in D$, we have that $z$ is non-negative for a.e. $(x,t)$. Moreover, we have:

[TABLE]

so that $z$ is essentially bounded. Fix $\varepsilon>0$, and take as $\psi(x,t)$ the solution of the following Neumann problem

[TABLE]

where the $\varepsilon$ is added to ensure the conjugate diffusion coefficient is strictly positive, and $F$ is an arbitrary smooth bounded function in $D$. Note that (3.19) is the conjugate parabolic equation. From [35], there exists a unique solution $\psi^{\varepsilon}\in W_{2}^{2,1}(D)$ of the problem (3.19)-(3.21). We will use the arbitrariness of $F$ to obtain that $B-\tilde{B}=0$ a.e. $(x,t)\in D$. Note that through (3.19), we can rewrite (3.18):

[TABLE]

Thus our goal will be attained if we have an energy estimate on $\mathcal{L}\psi$ for solutions of (3.19). For simplicity, we obtain energy estimates through the second order parabolic equation, which will give analogous estimates for the conjugate parabolic equation by reversing the time variable. Let $z^{\varepsilon}(x,t)=z(x,t)+\varepsilon$, and for simplicity we will omit the superscript. Multiply the non-conjugate version of (3.19) by $\psi_{xx}$ and integrate it over $D_{t}:=(0,\ell)\times(0,t)$ to get

[TABLE]

Due to (3.20), the second integral on the right hand side disappears. We can transform various terms on the right hand sign as follows:

[TABLE]

Using the above, and returning to (3.23), we get that:

[TABLE]

We now estimate the terms on the right hand side using Cauchy inequality with $\epsilon>0$ and properties of given functions, and absorbing terms to the left hand side, we have:

[TABLE]

From Theorem 2.3. Chapter 1 of [35] it follows that $\psi$ will have a uniform bound in $L_{\infty}(D)$, which we will denote as $\bar{C}$. Letting now $y(t)=\int\limits_{0}^{t}\int\limits_{0}^{\ell}\psi_{x}^{2}\,dx\,d\tau$ from (3.25) we deduce

[TABLE]

where

[TABLE]

By Gronwall’s Inequality (e.g. Lemma 5.5, Chapter 2, [35]), we deduce from the above differential inequality that

[TABLE]

so that by (3.25),

[TABLE]

The first of the above inequalities implies that

[TABLE]

Now, since $\psi_{t}=az\psi_{xx}+z(a_{x}-b)\psi_{x}+cz\psi+F$, we have

[TABLE]

where the constants $C_{3}$ and $C_{4}$ depend on $\bar{b},a_{0},T,\ell$ and $L_{\infty}$-norms of $a,a_{x},b,c$. Combining all the estimations we have the following desired energy estimate for $\psi\in W_{2}^{2,1}(D)$:

[TABLE]

where the constants $C_{5}$ and $C_{6}$ are independent of $\epsilon$, and depend on $\bar{b},a_{0},T,\ell$ and $L_{\infty}$-norms of $a,a_{x},b,c$.

For the rest of the proof, any constant depending on the bounded data, domain, $F$, or the uniform bound on $\psi$ will be referred to as $C^{*}$.We can now observe that

[TABLE]

as $\varepsilon\to 0$. Therefore, (3.22) now implies

[TABLE]

Since $F$ is arbitrary, the above equality gives that $B(x,t,v(x,t))=\tilde{B}(x,t,\tilde{v}(x,t))$ a.e. $(x,t)\in D$. This implies $\beta(v(x,t))=\beta(\tilde{v}(x,t))$, a.e. (x,t)\in D\leavevmode\nobreak\s.t.$\leavevmode\nobreak\ v(x,t)\neq v^{j},j=1,...,m$. Due to the fact that $\beta$ is strictly increasing, we have $v(x,t)=\tilde{v}(x,t)$ a.e. $(x,t)$. Thus $v$ and $\tilde{v}$ are the same solution to (3.17) $\square$

Corollary 1**.**

If $v$ is weak solution, the sets $\{(x,t)\in D|v=v^{j}\},j=1,...,m$ have 2-dimensional measure 0.

Indeed, due to uniqueness of the weak solution, for any two representatives $B_{1},B_{2}$ of the class $\mathscr{B}$ we have

[TABLE]

and by the Definition 1.1 this will be a contradiction if any of the $v^{j}$-level sets of the weak solution would have a positive 2-dimensional measure.

4. Proof of Main Results

4.1. $L_{\infty}(D)$ estimate for the discrete multiphase free boundary problem

In this section we prove $L_{\infty}(D)$ bound for the discrete PDE problem under the following reduced assumptions:

[TABLE]

and $a$ satisfies (1.7).

Theorem 4.1**.**

For $[g]_{n}\in\mathscr{G}_{R}^{n}$ and $n,m$ large enough, the discrete state vector $[v([g]_{n})]_{n}$ satisfies the following estimate:

[TABLE]

where $C_{\infty}$ is a constant independent of $n$ and $m$.

*Proof. * Fix $n$ arbitrarily large. Note $\max|v_{i}(0)|\leq\|\Phi\|_{L_{\infty}(0,\ell)}$. Consider a positive function $\gamma(x)\in C^{2}[0,\ell]$ satisfying

[TABLE]

Define $\gamma_{i}=\gamma(x_{i}),\leavevmode\nobreak\ i=\overline{0,m}$, and denote as $x^{i}$ the value in $[x_{i},x_{i+1}]$ that satisfies (by mean value theorem (MVT)) $\gamma(x_{i+1})-\gamma(x_{i})=\gamma^{\prime}(x^{i})h$. Transform the discrete state vector as

[TABLE]

System (3.4) can be rewritten as:

[TABLE]

We note

[TABLE]

Thus $w_{i}(0)=\gamma_{i}\Phi_{i},\quad i=\overline{0,m}$, and for $k=\overline{1,n}$,

[TABLE]

Furthermore, transform $w_{i}(k)$ as:

[TABLE]

where $\lambda$ satisfies

[TABLE]

and if $t^{k}\in[t_{k-1},t_{k}]$ satisfies through the MVT that $e^{\lambda t_{k}}-e^{\lambda t_{k-1}}=\lambda e^{\lambda t^{k}}\tau$, then

[TABLE]

So $u_{i}(0)=w_{i}(0)=\gamma_{i}\Phi_{i},\quad i=\overline{0,m}$, and for $k=\overline{1,n}$, the vector $u(k)$ satisfies the system

[TABLE]

Now fix $k_{1}\leq n$, and define the following sets of indexes for convenience:

[TABLE]

Unless confusion may arise, we omit the subscript to $\mathscr{M}_{k_{1}}$. It is clear that

[TABLE]

If $u_{i}(k)\leq 0$ in $\mathscr{M}$, then $\max\limits_{\mathscr{M}}u_{i}(k)\leq 0$. Suppose that there exists $(i,k)$ such that $u_{i}(k)>0$. Then $\max\limits_{\mathscr{M}}u_{i}(k)>0$. Let $(i^{*},k^{*})\in\mathscr{M}$ be such that $u_{i^{*}}(k^{*})=\max\limits_{\mathscr{M}}u_{i}(k)$.

If $(i^{*},k^{*})\in\mathscr{X}_{0}$, then $u_{i^{*}}(k^{*})=\max\limits_{i}\gamma_{i}\Phi_{i}\leq\max\limits_{i}\Phi_{i}\leq\max\limits_{[0,\ell]}\Phi(x)$.

If $(i^{*},k^{*})\in\mathscr{T}_{m}$, then $i^{*}=m,\leavevmode\nobreak\ u_{m-1,x}(k^{*})\geq 0$ and we can choose $h$ small enough that $\gamma_{m-1,x}=\gamma^{\prime}(x^{m-1})\in(-\frac{3}{2},-\frac{1}{2})$ so that

[TABLE]

We can see that:

[TABLE]

for sufficiently small h. Thus we have:

[TABLE]

If $(i^{*},k^{*})\in\mathscr{T}_{0}$, then $i^{*}=0,u_{0\bar{t}}(k^{*})\geq 0$, $u_{0x}(k^{*})\leq 0$. Notice that $\left(\frac{1}{\gamma_{0}}\right)_{x}=-\frac{1}{\gamma_{0}\gamma_{1}}\gamma_{0x}$. Note $\gamma_{0x}=\gamma^{\prime}(x^{0})$, so for $h$ small enough, we can ascertain $\gamma_{0x}=\gamma^{\prime}(x^{0})\in(\frac{4(\|b\|_{L_{\infty}(D)}+\|c\|_{L_{\infty}(D)})}{a_{0}}+\frac{1}{2},\frac{4(\|b\|_{L_{\infty}(D)}+\|c\|_{L_{\infty}(D)})}{a_{0}}+\frac{3}{2})$. It follows

[TABLE]

Since the third term in the parenthesis on the left hand side is positive, we only consider first two terms. We can see that:

[TABLE]

Thus we have:

[TABLE]

If $(i^{*},k^{*})\in\mathscr{N}$, then $u_{i^{*}\bar{t}}(k^{*})\geq 0$, \leavevmode\nobreak\ u_{i^{*}x\bar{x}}(k^{*})=\frac{1}{h^{2}}\big{(}u_{i^{*}+1}(k^{*})-2u_{i^{*}}(k^{*})+u_{i^{*}-1}(k^{*})\big{)}\leq 0. For $(i,k)\in\mathscr{N}$, the corresponding equation in (4.7) is equivalent to

[TABLE]

Define the sets

[TABLE]

And it’s clear $\mathscr{N}=\mathscr{N}_{+}\cup\mathscr{N}_{-}$. Suppose $(i^{*},k^{*})\in\mathscr{N}_{+}$. Then owing to (4.7) since $u_{i^{*}x}(k^{*})\leq 0$, for sufficiently small h we can write

[TABLE]

If instead $(i^{*},k^{*})\in\mathscr{N}_{-}$, then we can use (4.8), the fact that $\leavevmode\nobreak\ u_{i^{*}\bar{x}}(k^{*})\geq 0$ and that

[TABLE]

for sufficiently small h, where $\tilde{x}$ is the value that satisfies the mean value theorem to achieve again (4.9). Therefore, (4.9) is achieved in any case. We can choose $\tau$ so small that $e^{-\lambda(t_{k}-t^{k})}>\frac{1}{2},\quad\forall k$. The coefficient in front of $u_{i^{*}}(k^{*})$ in (4.9) can be estimated as follows:

[TABLE]

due to definitions of $\lambda$ and $\gamma(x)$. Then by (4.1), it is the case that the coefficient of $u_{i^{*}}(k^{*})$ is positive independently of $i^{*},k^{*}$. Therefore,

[TABLE]

We can put together the obtained estimations to deduce that for $(i,k)\in\mathscr{M}_{k_{1}}$,

[TABLE]

with $A=\max\{1,\leavevmode\nobreak\ 4a_{0}^{-1},\leavevmode\nobreak\ 2\bar{b}^{-1}\}$. But because $u_{i}(k)=\gamma_{i}e^{-\lambda t_{k}}v_{i}(k)$, we have the following uniform upper bound for the discrete state vector:

[TABLE]

for $\forall(i,k)\in\mathscr{M}_{k_{1}}$. In a fully analogous manner, we arrive at a uniform lower bound for the discrete state vector:

[TABLE]

for $\forall(i,k)\in\mathscr{M}_{k_{1}}$. Combining the uniform upper and lower bounds imply (4.1) up to $k_{1}$. But $k_{1}$ was arbitrary in $1,\ldots,n$. Theorem is proved. $\square$

4.2. $W_{2}^{1,1}(D)$ estimate for the discrete multiphase free boundary problem

Theorem 4.2**.**

For $[g]_{n}\in\mathscr{G}_{R}^{n}$ and $n,m$ large enough, the discrete state vector $[v([g]_{n})]_{n}$ satisfies the following estimate:

[TABLE]

where $\tilde{C}_{\infty}$ is a constant independent of $n$ and $m$.

*Proof. *Consider $n$ and $m$ large enough that Theorem 4.1 is satisfied. In (1.17), choose $\eta=2\tau v_{\bar{t}}(k)$. Using (3.13), write $(b_{n}(v_{i}(k)))_{\bar{t}}=\zeta_{k}^{i}v_{i\bar{t}}(k)$. Also, use the fact that

[TABLE]

Using the above equality, and the lower bound for $a(x,t)$, we thus have

[TABLE]

We will now look to estimate the three summation terms on the right hand side of (4.11). By using summation by parts and Cauchy inequality with $\epsilon>0$ we get:

[TABLE]

We will also use the fact that

[TABLE]

Estimating the other two summation terms on the right-hand side of (4.11) via Cauchy Inequality with $\varepsilon>0$ and by recalling (1.16), we have:

[TABLE]

By absorbing several terms on the right hand side of (4.12) to the left-hand side, and further bounding the right hand side we get

[TABLE]

$\forall k=\overline{1,n}$. Perform summation of (4.13) for $k$ from $1$ to $q,\leavevmode\nobreak\ \leavevmode\nobreak\ 2\leq q\leq n$. The second and third terms on the left-hand side telescope, and we obtain:

[TABLE]

We can estimate the right hand side further and use (1.9) to get

[TABLE]

Similarly,

[TABLE]

We also have:

[TABLE]

Since

[TABLE]

from (4.15) it follows that

[TABLE]

Similarly, we have:

[TABLE]

From [5], in a similar fashion to above, we have:

[TABLE]

where $C$ is a constant independent of $n$. Using Cauchy Inequality with $\varepsilon=\frac{\bar{b}}{8}$, we have:

[TABLE]

Use the summation by parts technique on the $p$ and $g$ sums:

[TABLE]

In view of (4.16) and the above estimates, (4.14) yields

[TABLE]

Applying Cauchy-Schwartz inequality to Steklov averages, by emloying Morrey’s inequality ([35]), for sufficiently small h we have the following estimations:

[TABLE]

with last two estimations being extended to similar terms

[TABLE]

Applying the results in (4.18), along with $L_{\infty}$-estimate (4.1) to (4.17), and taking into account that $q=\overline{1,n}$ is arbitrary we derive

[TABLE]

where $\tilde{C}_{\infty}$ is a constant independent of $n,m$. If

[TABLE]

then we absorb the extra term to the left-hand side and (4.10) follows. If not, we partition $[0,T]$ into finitely many intervals such that in each interval $I$, (4.20) is satisfied with integral along $I$, and $T$ in the second term replaced with the interval length $|I|$. Therefore, energy estimate (4.10) holds in each interval segment, and by adding finitely many inequalities, (4.10) in the whole segment follows. $\square$

4.3. Existence of Optimal Control and Convergence of Discrete Control Problems

Theorem 4.3**.**

Let $\{[g]_{n}\}$ be a sequence in $\mathscr{G}_{R}^{n}$ such that the sequence of interpolations $\{\mathscr{P}_{n}([g]_{n})\}$ converges weakly to $g\in W_{2}^{1}[0,T]$. Then the whole sequence of interpolations $\{\hat{v}^{\tau}\}$ of the associated discrete state vectors converges weakly in $W_{2}^{1,1}(D)$ to the unique weak solution $v=v(x,t;g)$ of the multiphase parabolic free boundary problem (1.2)-(1.5).

*Proof. * Having estimates (4.1),(4.2), the proof is pursued similar to the proof of Theorem 5 in [5]. By the definitions of the interpolations in (1.19) we have that there is a subsequence of $\{\hat{v}^{\tau}\}$ that converges weakly in $W_{2}^{1,1}(D)$ to some function $v\in W_{2}^{1,1}(D)\cap L_{\infty}(D)$, strongly in $L_{2}(D)$, and a further subsequence that converges to $v$ pointwise almost everywhere in $D$. We also have equivalence of $\{v^{\tau}\}$ and $\{\hat{v}^{\tau}\}$ in $W_{2}^{1,0}(D)$, and equivalence of $\{v^{\tau}\}$ and $\{\tilde{v}\}$ in $L_{2}(D)$. Accordingly, $v^{\tau}\rightarrow v$ weakly in $W_{2}^{1,0}(D)$ and $\tilde{v}\rightarrow v$ strongly in $L_{2}(D)$ and pointwise a.e. on $D$ along a subsequence. Fix arbitrary $\psi\in W_{2}^{1,1}(D)$ with $\psi|_{t=T}=0$. Due to density of $C^{1}(\bar{D})$ in $W_{2}^{1,1}(D)$, without loss of generality we can consider $\psi\in C^{1}(\overline{D})$ and $\psi|_{t=T}=0$. Define $\psi_{i}(k)=\psi(x_{i},t_{k})$, and consider the interpolations:

[TABLE]

It is readily checked that $\psi^{\tau},\psi_{x}^{\tau},\psi_{t}^{\tau}$ converge uniformly on $\overline{D}$ as $n,m\rightarrow\infty$ to the functions $\psi,\psi_{x},\psi_{t}$ respectively. For each $k$ in (1.17) as satisfied by the discrete state vector $[v([g]_{n})]_{n}$, choose $\eta_{i}=\tau\psi_{i}(k),\leavevmode\nobreak\ i=0,...,m$ and sum all equalities (1.17) over $k=1,\ldots,n$. The resulting expression is as follows:

[TABLE]

The first term is transformed through summation by parts as in [5], and using the interpolations, (4.22) becomes the following integral identity:

[TABLE]

From [5], we have that $b_{n}(\tilde{v})$ has a subsequence such that both $b_{n}(\tilde{v})$ and $b_{n}(\tilde{\Phi})$ converge weakly in $L_{2}(D)$ and $L_{2}[0,\ell]$, respectfully, to functions of type $\mathscr{B}$, which we will denote as $\tilde{b}(x,t)$ and $\tilde{b}_{0}(x)$. We also have that the last integral tends to [math] due to absolute continuity of the integral. Using the convergence properties of the interpolations, due to weak convergence of $\{\hat{v}^{\tau}\}$, equivalence of $\{\hat{v}^{\tau}\}$ and $\{v^{\tau}\}$, and uniform convergence of $\{\psi^{\tau}_{x}\}$, passing to the limit as $n\to+\infty$ we get:

[TABLE]

Since $\tilde{b}(x,t)$ and $\tilde{b}_{0}(x)$ are both of type $\mathscr{B}$, and by use of Mazur’s lemma, we deduce as in [5] that $\tilde{b}(x,t)=B(x,t,v(x,t))$ $\tilde{b}_{0}(x)=B(x,0,\Phi(x))$ a.e on $D$ and $(0,\ell)$ respectfully. This implies that $v$ is a weak solution in the sense of Definition 1.2. By Lemma 3.4, this implies $v$ is the only solution of the problem, and hence the only limit point of the sequence $\{\hat{v}^{\tau}\}$. $\square$

Having estimates (4.1),(4.2) and compactness Theorem 4.3, the completion of the proof of Theorems 2.1 and 2.2 coincides with the proof given in [5], through compactness arguments and proving weak continuity of cost functional $\mathscr{J}$, and verification of the conditions of Theorem 3.2 (see Lemmas A, B and C in [5]).