Recovering the initial condition of parabolic equations from lateral Cauchy data via the quasi-reversibility method

ABSTRACT

We consider the problem of computing the initial condition for a general parabolic equation from the Cauchy lateral data. The stability of this problem is well-known to be logarithmic. In this paper, we introduce an approximate model, as a coupled linear system of elliptic partial differential equations. Solution to this model is the vector of Fourier coefficients of the solutions to the parabolic equation above. This approximate model is solved by the quasi-reversibility method. We will prove the convergence for the quasi-reversibility method as the measurement noise tends to 0. The convergent rate is Lipschitz. We present the implementation of our algorithm in details and verify our method by showing some numerical examples.

KEYWORDS:

• Initial condition
• parabolic equation
• truncated Fourier series
• approximation
• convergence
• quasi-reversibility method

2010 MATHEMATICS SUBJECT CLASSIFICATIONS:

• 35R30
• 35K15

1. Introduction

Let $d\ge 2$ be the spatial dimension and Ω be a open and bounded domain in ${\mathbb{R}}^d$. Assume that $\mathrm{\partial }\mathrm{\Omega }$ is smooth. Let (1) $A=(a_{ij}{)}_{i,j=1}^d\in C^2({\mathbb{R}}^d,{\mathbb{R}}^{d\times d})\cap L^{\mathrm{\infty }}({\mathbb{R}}^d,{\mathbb{R}}^{d\times d})$(1) satisfy the following conditions

1. A is symmetric; i.e. $A^T(\mathbf{x})=A(\mathbf{x})$ for all $\mathbf{x}\in {\mathbb{R}}^d;$

2. A is uniformly elliptic; i.e. there exists a positive number μ such that (2) $A(\mathbf{x})\xi \cdot \xi \ge \mu |\xi {|}^2\mathrm{for}\mathrm{all}\mathbf{x},\xi =({\xi }_1,\dots ,{\xi }_d)\in {\mathbb{R}}^d.$(2)

Let $\mathbf{b}=(b_1,b_2,\dots ,b_d)\in C^1({\mathbb{R}}^d,{\mathbb{R}}^d)\cap L^{\mathrm{\infty }}({\mathbb{R}}^d,{\mathbb{R}}^d)$ and $c\in C^1({\mathbb{R}}^d,\mathbb{R})\cap L^{\mathrm{\infty }}({\mathbb{R}}^d,\mathbb{R})$. Define (3) $\mathcal{L}v=\mathrm{Div}(A\mathrm{\nabla }v)+\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }v(\mathbf{x})+c(\mathbf{x})v$(3) for all functions $v\in C^2({\mathbb{R}}^d)$. Consider the initial value problem (4) $\begin{array}{rl}{u}_t(\mathbf{x},t)& =\mathcal{L}u(\mathbf{x},t),\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ u(\mathbf{x},0)& =f(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,\end{array}$(4) where $f\in L^2({\mathbb{R}}^d)$ represents an initial source with support compactly contained in Ω. We refer the reader to the books [1,2]. The main aim of this paper is to solve the following problem.

problem 1.1

Let T>0. Given the Cauchy boundary data (5) $F(\mathbf{x},t)=u(\mathbf{x},t)\mathrm{and}G(\mathbf{x},t)=A\mathrm{\nabla }u(\mathbf{x},t)\cdot \nu$(5) for $\mathbf{x}\in \mathrm{\partial }\mathrm{\Omega },t\in [0,T],$ determine the function $f(\mathbf{x}),$ $\mathbf{x}\in \mathrm{\Omega }$. Here ν is the outward normal to $\mathrm{\partial }\mathrm{\Omega }.$

Problem 1.1 is the problem of recovering the initial condition of the parabolic equation from the lateral Cauchy data. This problem has many real-world applications; e.g. determine the spatially distributed temperature inside a solid from the boundary measurement of the heat and heat flux in the time domain [3]; identify the pollution on the surface of the rivers or lakes [4]; effectively monitor the heat conductive processes in steel industries, glass and polymer forming and nuclear power station [5]. Due to its realistic applications, this problem has been studied intensively. The uniqueness of Problem 1.1 is well-known, see [6]. Also, it can be reduced from the logarithmic stability results in [3,5]. The natural approach to solve this problem is the optimal control method; that means, minimize a mismatch functional. The proof of the convergence of the optimal control method to the true solution to these inverse problems is challenging and is omitted. One of our contributions to the field is the convergence of the quasi-reversibility method, which our method is relied on, as the measurement noise tends to 0.

Related to the inverse problem in the current paper, the problem of recovering the initial conditions for hyperbolic equation is very interesting since it arises in many real-world applications. For instance the problems thermo and photo acoustic tomography play the key roles in bio-medical imaging. We refer the reader to some important works in this field [7–9]. Applying the Fourier transform, one can reduce the problem of reconstructing the initial conditions for hyperbolic equations to some inverse source problems for the Helmholtz equation, see [10–14] for some recent results.

In this paper, we employ the technique developed by our own research group. The main point of this technique is to derive an approximate model for the Fourier coefficients of the solution to the governing partial differential equation. This technique was first introduced in [15]. This approximate model is a system of elliptic equations. It, together with Cauchy boundary data, is solved by the quasi-reversibility method. This approach was used to solve an inverse source problem for Helmholtz equation [10] and to inverse the Radon transform with incomplete data [16]. Especially, Klibanov et al. [17] used the convexification method, a stronger version of this technique, to compute numerical solutions to the nonlinear problem of electrical impedance tomography with restricted Dirichlet-to-Neumann map data. It is remarkable mentioning that the numerical solutions in [17] due to the convexification method are impressive.

As mentioned in the previous paragraph, we employ the quasi-reversibility method to solve an approximate model for Fourier coefficients of the solution to (4(4) $\begin{array}{rl}{u}_t(\mathbf{x},t)& =\mathcal{L}u(\mathbf{x},t),\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ u(\mathbf{x},0)& =f(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,\end{array}$(4) ). This method was first introduced by Lattés and Lions [18]. It is used to computed numerical solutions to ill-posed problems for partial differential equations. Due to its strength, since then, the quasi-reversibility method attracts the great attention of the scientific community see e.g. [19–28]. We refer the reader to [29] for a survey on this method. The solution of the approximate model in the previous paragraph due to the quasi-reversibility method is called regularized solution in the theory of ill-posed problems [30]. A question arises immediately about the convergence of the quasi-reversibility method: whether or not the regularized solutions obtained by the quasi-reversibility method converges to the true solution of our system of partial differential equations as the noise tends to 0. The affirmative answer to this question is obtained using a general Carleman estimate. Moreover, we employ a Carleman estimate to prove that the convergence rate is Lipschitz. It is important mentioning that in the celebrate paper [31], Bukhgeim and Klibanov discovered the use of Carleman estimate in studying inverse problems for all three main types of partial differential equations.

The paper is organized as follows. In Section 2, we describe our approach and propose an algorithm to solve Problem 1.1. In Section 3, we employ prove a Carleman estimate. Then, in Section 4, we study the convergence of the quasi-reversibility method as the noise tends to 0. Finally, in Section 5, we present all details about the numerical implementation and then show some numerical results from highly noisy simulated data.

2. The algorithm to solve Problem 1.1

We will employ the following basis to introduce an approximation model.

2.1. An orthonormal basis of $L^2(0,T)$ and the truncated Fourier series

For each n>1, define a complete sequence $\{{\varphi }_n{\}}_{n=1}^{\mathrm{\infty }}$ in $L^2(0,T)$ with (6) ${\varphi }_n(t)=(t-t_0{)}^{n-1}\exp (t-t_0),$(6) where $t_0=T/2.$ Using the Gram–Schmidt orthonormalization for the sequence $\{{\varphi }_n{\}}_{n=1}^{\mathrm{\infty }}$, we can construct an orthonormal basis of $L^2(0,T),$ named as $\{{\mathrm{\Psi }}_n{\}}_{n=1}^{\mathrm{\infty }}$. For each n, the function ${\mathrm{\Psi }}_n(t)$ takes the form (7) ${\mathrm{\Psi }}_n(t)=P_{n-1}(t-t_0)\exp (t-t_0),$(7) where $P_{n-1}$ is a polynomial of the $(n-1)\mathrm{th}$ order. For each $\mathbf{x}\in \mathrm{\Omega },$ we consider $u(\mathbf{x},\cdot )$ as a function with respect to t. The Fourier series of this function is (8) $u(\mathbf{x},t)=\sum _{n=1}^{\mathrm{\infty }}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),$(8) where (9) $u_n(\mathbf{x})={\int }_0^{T}u(\mathbf{x},t){\mathrm{\Psi }}_n(t)\mathrm{d}t,n=1,2,\dots .$(9) Fix a positive integer N. We truncate the Fourier series in (8(8) $u(\mathbf{x},t)=\sum _{n=1}^{\mathrm{\infty }}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),$(8) ). The function $u(\mathbf{x},t)$ is approximated by (10) $u(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),\mathbf{x}\in \mathrm{\Omega },t\in [0,T].$(10) In this context, the partial derivative with respect to t of $u(\mathbf{x},t)$ is approximated by (11) $u_t(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n^{\prime }(t)$(11) for all $\mathbf{x}\in \mathrm{\Omega }$ and $t\in (0,T).$

To reconstruct the wave field $u(\mathbf{x},t)$, we compute the Fourier coefficients $u_n(\mathbf{x})$, $1\le n\le N$. It is obvious that (10(10) $u(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),\mathbf{x}\in \mathrm{\Omega },t\in [0,T].$(10) ) and (11(11) $u_t(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n^{\prime }(t)$(11) ) play crucial roles in this step. We; therefore, require that the function ${\mathrm{\Psi }}_n^{\prime }$ cannot be identically 0. The usual ‘sin and cosine’ basis of the Fourier transform does not meet this requirement while it is not hard to verify from (7(7) ${\mathrm{\Psi }}_n(t)=P_{n-1}(t-t_0)\exp (t-t_0),$(7) ) that the basis $\{{\mathrm{\Psi }}_n{\}}_{n=1}^{\mathrm{\infty }}$ does. The basis $\{{\mathrm{\Psi }}_n{\}}_{n=1}^{\mathrm{\infty }}$ was first introduced in [15]. Then, this basis was successfully used to solve several important inverse problems, including the inverse source problem for Helmholtz equations [10], inverse X-ray tomographic problem in incomplete data [16] and the nonlinear inverse problem of electrical impedance tomography with restricted Dirichlet to Neumann map data, see [17].

2.2. An approximate model

We introduce in this subsection a coupled system of elliptic partial differential equations without the presence of the unknown function $f(\mathbf{x})$. Plugging (10(10) $u(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),\mathbf{x}\in \mathrm{\Omega },t\in [0,T].$(10) ) and (11(11) $u_t(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n^{\prime }(t)$(11) ) into (4(4) $\begin{array}{rl}{u}_t(\mathbf{x},t)& =\mathcal{L}u(\mathbf{x},t),\mathbf{x}\in {\mathbb{R}}^d,t>0,\\ u(\mathbf{x},0)& =f(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,\end{array}$(4) ), we have (12) $\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n^{\prime }(t)=\sum _{n=1}^N\mathcal{L}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t)$(12) for all $\mathbf{x}\in \mathrm{\Omega }$ and $t\in [0,T].$ For each $m\in \{1,\dots ,N\}$, multiply ${\mathrm{\Psi }}_m(t)$ to both sides of (12(12) $\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n^{\prime }(t)=\sum _{n=1}^N\mathcal{L}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t)$(12) ) and then integrating the obtained equation with respect to t, we obtain (13) $\sum _{n=1}^N({\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n^{\prime }(t)\mathrm{d}t)u_n(\mathbf{x})=\sum _{n=1}^N({\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n(t)\mathrm{d}t)\mathcal{L}{u}_n(\mathbf{x})$(13) for all $\mathbf{x}$ in Ω. Denote by (14) $s_{mn}={\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n^{\prime }(t)\mathrm{d}t$(14) and note that (15) ${\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n(t)\mathrm{d}t=\{\begin{array}{ll}1& \mathrm{if}m=n,\\ 0& \mathrm{if}m\ne n.\end{array}$(15) We rewrite (13(13) $\sum _{n=1}^N({\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n^{\prime }(t)\mathrm{d}t)u_n(\mathbf{x})=\sum _{n=1}^N({\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n(t)\mathrm{d}t)\mathcal{L}{u}_n(\mathbf{x})$(13) ) as (16) $\sum _{n=1}^N{s}_{mn}{u}_n(\mathbf{x})=\mathcal{L}{u}_m(\mathbf{x}),\mathbf{x}\in \mathrm{\Omega },m=1,2,\dots ,N.$(16) Denote (17) $U(\mathbf{x})=(u_1(\mathbf{x}),u_2(\mathbf{x}),\dots ,u_N(\mathbf{x}){)}^T.$(17) It follows from (16(16) $\sum _{n=1}^N{s}_{mn}{u}_n(\mathbf{x})=\mathcal{L}{u}_m(\mathbf{x}),\mathbf{x}\in \mathrm{\Omega },m=1,2,\dots ,N.$(16) ) that (18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) where S is the $N\times N$ matrix whose $mn\mathrm{th}$ entry is given in (14(14) $s_{mn}={\int }_0^T{\mathrm{\Psi }}_m(t){\mathrm{\Psi }}_n^{\prime }(t)\mathrm{d}t$(14) ), $1\le m,n\le N$. Here, the operator $\mathcal{L}$ acting on the vector $U(\mathbf{x})$ is understood in the same manner as it acts on scalar valued function, see (3(3) $\mathcal{L}v=\mathrm{Div}(A\mathrm{\nabla }v)+\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }v(\mathbf{x})+c(\mathbf{x})v$(3) ).

On the other hand, due to (9(9) $u_n(\mathbf{x})={\int }_0^{T}u(\mathbf{x},t){\mathrm{\Psi }}_n(t)\mathrm{d}t,n=1,2,\dots .$(9) ) and (5(5) $F(\mathbf{x},t)=u(\mathbf{x},t)\mathrm{and}G(\mathbf{x},t)=A\mathrm{\nabla }u(\mathbf{x},t)\cdot \nu$(5) ), the vector $U(\mathbf{x})$ satisfies the boundary conditions (19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) (20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) for all $\mathbf{x}\in \mathrm{\partial }\mathrm{\Omega }.$

Remark 2.1

From now on, we consider $\tilde{F}$ and $\tilde{G}$ as our ‘indirect’ boundary data. This is acceptable since these two functions can be computed directly by the algebraic formulas (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ).

Finding a vector $U(\mathbf{x})$ satisfying Equation (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ) and constraints (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) is the main point in our numerical method to find the function $f(\mathbf{x})$. In fact, having $U(\mathbf{x})=(u_1(\mathbf{x}),\dots ,u_2(\mathbf{x}),\dots ,u_N(\mathbf{x}))$ in hand, we can compute the function $u(\mathbf{x},t)$ via (10(10) $u(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),\mathbf{x}\in \mathrm{\Omega },t\in [0,T].$(10) ). The desired function $f(\mathbf{x})$ is given by $u(\mathbf{x},t=0).$

Due to the truncation step in (10(10) $u(\mathbf{x},t)=\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),\mathbf{x}\in \mathrm{\Omega },t\in [0,T].$(10) ), Equation (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ) is not exact. We call it an approximate model. Solving it, together with the ‘over-determined’ boundary conditions (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ), for the Fourier coefficients $(u_n(\mathbf{x}){)}_{n=1}^N$ of $u(\mathbf{x},t)$, $\mathbf{x}\in \mathrm{\Omega }$, $t\in [0,T]$, might not be rigorous. In fact, proving the ‘accuracy’ of (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ) when $N\to \mathrm{\infty }$ is extremely challenging and is out of the scope of this paper. However, we experience in many earlier works that the solution of (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) well approximates Fourier coefficients of the function $u(\mathbf{x},t)$, leading to good solutions of variety kinds of inverse problems, see [10,16,17,32].

Remark 2.2

The choice of N

On $\mathrm{\Omega }=(-2,2{)}^2$, we arrange $81\times 81$ grid points $\{(x_i,y_j):1\le i,j\le 81\}$. In Figure 1 displays the functions of $u(\mathbf{x},t)$ and its approximation $\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t)$ where $u(\mathbf{x},t)$ is the true solution of the forward problem and $u_n(\mathbf{x})$, $n=1,\dots ,N$, is computed using (9(9) $u_n(\mathbf{x})={\int }_0^{T}u(\mathbf{x},t){\mathrm{\Psi }}_n(t)\mathrm{d}t,n=1,2,\dots .$(9) ). This numerical experiment suggests us to take N=30. It is worth mentioning that when $N\le 25$, the numerical solutions are not satisfactory, when N=30, numerical results are quite accurate regardless the high noise levels and when $N\ge 35,$ the computation is time-consuming.

Figure 1. The function $u(\mathbf{x},t=0)$ (dash-dot) and its approximation $\sum _{n=1}^N{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t=0)$ (solid) at the points numbered from 900 to 1050. These functions are taken from Test 4 in Section 5.2. It is evident that the larger N, the better approximation for the function u is obtained by the $N\mathrm{th}$ partial sum of the Fourier series in (8(8) $u(\mathbf{x},t)=\sum _{n=1}^{\mathrm{\infty }}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),$(8) ): (a) N=10, (b) N=20, (c) N=30.

2.3. The quasi-reversibility method

As mentioned, our method to solve Problem 1.1 is based on a numerical solver for (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ). We do so by employing the quasi-reversibility method; that means, we minimize the functional (21) $J_{ϵ}(U)={\int }_{\mathrm{\Omega }}|\mathcal{L}U(\mathbf{x})-SU(\mathbf{x}){|}^2\mathrm{d}\mathbf{x}+ϵ\Vert U{\Vert }_{H^3(\mathrm{\Omega })}^2$(21) subject to the constraints (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ). Here ε is a positive number serving as a regularization parameter. Impose the condition that the set of admissible data (22) $H=\{V\in H^3(\mathrm{\Omega }{)}^N:V{|}_{\mathrm{\partial }\mathrm{\Omega }\times [0,T]}=\tilde{F}\mathrm{and}A\mathrm{\nabla }V\cdot \nu {|}_{\mathrm{\partial }\mathrm{\Omega }\times [0,T]}=\tilde{G}\}$(22) is nonempty, where $\tilde{F}$ and $\tilde{G}$ are our indirect data, see Remark 2.1, defined in (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ). The result below guarantees the existence and uniqueness for the minimizer of $J_{ϵ}$, $ϵ>0$.

Proposition 2.3

Assume that the set of admissible data H, defined in (22(22) $H=\{V\in H^3(\mathrm{\Omega }{)}^N:V{|}_{\mathrm{\partial }\mathrm{\Omega }\times [0,T]}=\tilde{F}\mathrm{and}A\mathrm{\nabla }V\cdot \nu {|}_{\mathrm{\partial }\mathrm{\Omega }\times [0,T]}=\tilde{G}\}$(22) ), is nonempty. Then, for all $ϵ>0,$ the functional $J_{ϵ}$ admits a unique minimizer satisfying (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ). This minimizer is called the regularized solution to (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ).

Proof.

Proposition 2.3 is an analogue of [10, Theorem 3.1] whose proof is based on the Riesz representation theorem. An alternative method to prove this proposition is from the standard argument in convex analysis, see e.g. [28,33].

The minimizer of $J_{ϵ}$ in H is called the regularized solution of (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) obtained by the quasi-reversibility method.

The analysis above leads to Algorithm 1, which describes our numerical method to reconstruct the function $f(\mathbf{x})$, $\mathbf{x}\in \mathrm{\Omega }$. In the next section, we establish a new Carleman estimate. This estimate plays an important role in proving the convergence of the regularized solution, due to the quasi-reversibility method, to the true solution of (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) in Section 4 as the measurement noise and ε tend to 0.

3. A Carleman estimate for second order elliptic operators on general domains

Let the matrix A be as in (1(1) $A=(a_{ij}{)}_{i,j=1}^d\in C^2({\mathbb{R}}^d,{\mathbb{R}}^{d\times d})\cap L^{\mathrm{\infty }}({\mathbb{R}}^d,{\mathbb{R}}^{d\times d})$(1) ). The main aim of this section is to prove a Carleman estimate in a general domain Ω. Similar versions of Carleman estimate can be found in [17, Theorem 3.1] and [34, Lemma 5] when Ω is an annulus and [10, Theorem 4.1] and when Ω is a cube. In this paper, we will use the following estimate to derive the convergence of the quasi-reversibility method. It can be deduced from [6, Lemma 3, Chapter 4, § 1].

Without lost of generality, we can assume that (23) $\mathrm{\Omega }\subset \{\mathbf{x}=(x_1,x_2,\dots ,x_d):0<x_1+X^{-2}\sum _{i=2}^d{x}_i^2<1\}$(23) for some 0<X<1. Define the function (24) $\psi (\mathbf{x})=x_1+\frac{1}{2X^2}\sum _{i=1}^d{x}_i^2+\alpha ,0<\alpha <1/2.$(24) Using Lemma 3 in [6, Chapter 4, § 1] for the function $u\in C^2(\overline{\mathrm{\Omega }})$ that is independent of the time variable, we can find a constant ${\sigma }_0$ and a constant ${\sigma }_1$ (depending only on α and the entries $a_{ij}$, $1\le i,j\le d$, of the matrix A) such that for all $\lambda \ge {\sigma }_0$ and $p>{\sigma }_1$ (25) $\begin{array}{rl}\frac{\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2+{\lambda }^3{p}^4{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2& \le -\frac{C\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\\ & +C{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2+\mathrm{Div}U\end{array}$(25) for all $\mathbf{x}\in \mathrm{\Omega }$ where the vector U satisfies (26) $|U|\le C{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}(\frac{\lambda p}{X}|\mathrm{\nabla }u{|}^2+\frac{{\lambda }^3{p}^3}{X^3}{\psi }^{-2p-2}{u}^2).$(26) Applying (25(25) $\begin{array}{rl}\frac{\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2+{\lambda }^3{p}^4{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2& \le -\frac{C\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\\ & +C{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2+\mathrm{Div}U\end{array}$(25) ) and (26(26) $|U|\le C{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}(\frac{\lambda p}{X}|\mathrm{\nabla }u{|}^2+\frac{{\lambda }^3{p}^3}{X^3}{\psi }^{-2p-2}{u}^2).$(26) ), we have the lemma.

Lemma 3.1

Carleman estimate

Let $u\in C^2(\overline{\mathrm{\Omega }})$ satisfying (27) $u{|}_{\mathrm{\partial }\mathrm{\Omega }}=A\mathrm{\nabla }u\cdot \nu =0\mathit{\text{on}}\mathrm{\partial }\mathrm{\Omega },$(27) where ν the outward unit normal vector of $\mathrm{\partial }\mathrm{\Omega }.$ Then, there exist a positive number ${\sigma }_0$ and ${\sigma }_1,$ depending only on α and A, such that (28) $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}\end{array}$(28) for $\lambda >{\sigma }_0$ and $p>{\sigma }_1$. In particular, fixing $p>{\sigma }_1$, one can find $\lambda >{\sigma }_0$ such that (29) $\begin{array}{rl}& \lambda p{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le C{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(29)

Proof.

We claim that (30) $\mathrm{\nabla }u(\mathbf{x})=0\mathrm{for}\mathrm{all}\mathbf{x}\in \mathrm{\partial }\mathrm{\Omega }.$(30) In fact, assume that $\mathrm{\nabla }u(\mathbf{x})\ne 0$ at some points $\mathbf{x}\in \mathrm{\partial }\mathrm{\Omega }.$ Since $u(\mathbf{x})=0$ on $\mathrm{\partial }\mathrm{\Omega }$, see (27(27) $u{|}_{\mathrm{\partial }\mathrm{\Omega }}=A\mathrm{\nabla }u\cdot \nu =0\mathit{\text{on}}\mathrm{\partial }\mathrm{\Omega },$(27) ), $\mathrm{\nabla }u(\mathbf{x})\cdot \tau (\mathbf{x})=0$ where $\tau (\mathbf{x})$ is any tangent vector to $\mathrm{\partial }\mathrm{\Omega }$ at the point $\mathbf{x}$. Thus, $\mathrm{\nabla }u(\mathbf{x})$ is perpendicular to $\mathrm{\partial }\mathrm{\Omega }$ at $\mathbf{x}$. In other words, $\mathrm{\nabla }u(\mathbf{x})=\theta \nu (\mathbf{x})$ for some nonzero scalar θ. We have $0=A(\mathbf{x})\mathrm{\nabla }u(\mathbf{x})\cdot \nu (\mathbf{x})=\theta A(\mathbf{x})\nu (\mathbf{x})\cdot \nu (\mathbf{x})$, which is a contradiction to (2(2) $A(\mathbf{x})\xi \cdot \xi \ge \mu |\xi {|}^2\mathrm{for}\mathrm{all}\mathbf{x},\xi =({\xi }_1,\dots ,{\xi }_d)\in {\mathbb{R}}^d.$(2) ).

Integrating both sides of (25(25) $\begin{array}{rl}\frac{\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2+{\lambda }^3{p}^4{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2& \le -\frac{C\lambda p}{X^2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\\ & +C{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2+\mathrm{Div}U\end{array}$(25) ), we have (31) $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le -\frac{C\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\mathrm{d}\mathbf{x}+C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(31) Here, the term ${\displaystyle {\int }_{\mathrm{\Omega }}\mathrm{Div}U\mathrm{d}\mathbf{x}}$ is dropped because it vanishes due the divergence theorem, (27(27) $u{|}_{\mathrm{\partial }\mathrm{\Omega }}=A\mathrm{\nabla }u\cdot \nu =0\mathit{\text{on}}\mathrm{\partial }\mathrm{\Omega },$(27) ) and (30(30) $\mathrm{\nabla }u(\mathbf{x})=0\mathrm{for}\mathrm{all}\mathbf{x}\in \mathrm{\partial }\mathrm{\Omega }.$(30) ) Using the inequality $|ab|\le \lambda p{a}^2+(1/2\lambda p)b^2$ (32) $\begin{array}{rl}& -\frac{C\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\mathrm{d}\mathbf{x}\\ & \le \frac{C{\lambda }^2{p}^2}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}{u}^2\mathrm{d}\mathbf{x}+\frac{C}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(32) Combining (31(31) $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le -\frac{C\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\mathrm{d}\mathbf{x}+C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(31) ) and (32(32) $\begin{array}{rl}& -\frac{C\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}u\mathrm{Div}(A\mathrm{\nabla }u)\mathrm{d}\mathbf{x}\\ & \le \frac{C{\lambda }^2{p}^2}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}{u}^2\mathrm{d}\mathbf{x}+\frac{C}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(32) ), we obtain $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le \frac{C{\lambda }^2{p}^2}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}{u}^2\mathrm{d}\mathbf{x}+\frac{C}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}\\ & +C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$ Fixing $p>{\sigma }_1$ and choosing λ large such that the second term on the left-hand side dominates the first term on the right-hand side, we obtain $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$ The estimate (28(28) $\begin{array}{rl}& \frac{\lambda p}{X^2}{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{\psi }^{-2p-2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le C{\int }_{\mathrm{\Omega }}{\psi }^{p+2}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}\end{array}$(28) ) follows.

4. The convergence of the quasi-reversibility method

In this section, we continue to assume (23(23) $\mathrm{\Omega }\subset \{\mathbf{x}=(x_1,x_2,\dots ,x_d):0<x_1+X^{-2}\sum _{i=2}^d{x}_i^2<1\}$(23) ). Let ${\tilde{F}}^{\ast }$ and ${\tilde{G}}^{\ast }$ be the noiseless data for (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ), see Remark (2.1), respectively. The noisy data are denoted by ${\tilde{F}}^{\delta }$ and ${\tilde{G}}^{\delta }$. Here δ is the noise level. In this section, assume that there exists $\mathcal{E}\in H^3(\mathrm{\Omega }{)}^N$ such that

1. for all $\mathbf{x}\in \mathrm{\Omega }$ (33) $\mathcal{E}(\mathbf{x})={\tilde{F}}^{\delta }(\mathbf{x})-{\tilde{F}}^{\ast }(\mathbf{x})\mathrm{and}A(\mathbf{x})\mathrm{\nabla }\mathcal{E}(\mathbf{x})\cdot \nu (\mathbf{x})={\tilde{G}}^{\delta }(\mathbf{x})-{\tilde{G}}^{\ast }(\mathbf{x});$(33)

2. and the bound (34) $\Vert \mathcal{E}{\Vert }_{H^3(\mathrm{\Omega }{)}^N}<\delta \mathrm{as}\delta \to 0^{+}$(34) holds true.

The assumption about the existence of $\mathcal{E}$ satisfying (33(33) $\mathcal{E}(\mathbf{x})={\tilde{F}}^{\delta }(\mathbf{x})-{\tilde{F}}^{\ast }(\mathbf{x})\mathrm{and}A(\mathbf{x})\mathrm{\nabla }\mathcal{E}(\mathbf{x})\cdot \nu (\mathbf{x})={\tilde{G}}^{\delta }(\mathbf{x})-{\tilde{G}}^{\ast }(\mathbf{x});$(33) ) and (34(34) $\Vert \mathcal{E}{\Vert }_{H^3(\mathrm{\Omega }{)}^N}<\delta \mathrm{as}\delta \to 0^{+}$(34) ) is equivalent to the condition $inf\{\Vert \mathrm{\Phi }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}:\mathrm{\Phi }{|}_{\mathrm{\partial }\mathrm{\Omega }}={\tilde{F}}^{\delta }-{\tilde{F}}^{\ast },{\mathrm{\partial }}_{\nu }\mathrm{\Phi }{|}_{\mathrm{\partial }\mathrm{\Omega }}={\tilde{G}}^{\delta }-{\tilde{G}}^{\ast }\}<\delta .$ In this section, we establish the following result to study the accuracy of the quasi-reversibility method.

Theorem 4.1

Assume that $U^{\ast }$ is the function that satisfies (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ), (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) with $\tilde{F}$ and $\tilde{G}$ replaced by ${\tilde{F}}^{\ast }$ and ${\tilde{G}}^{\ast }$ respectively. Fix $ϵ>0.$ Let $U^{\delta }$ be the minimizer of $J_{ϵ}$ subject to constraints (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ) with $\tilde{F}$ and $\tilde{G}$ replaced by ${\tilde{F}}^{\delta }$ and ${\tilde{G}}^{\delta }$ respectively. Assume further that there is an ‘error’ function $\mathcal{E}$ in $H^3(\mathrm{\Omega }{)}^N$ satisfying (33(33) $\mathcal{E}(\mathbf{x})={\tilde{F}}^{\delta }(\mathbf{x})-{\tilde{F}}^{\ast }(\mathbf{x})\mathrm{and}A(\mathbf{x})\mathrm{\nabla }\mathcal{E}(\mathbf{x})\cdot \nu (\mathbf{x})={\tilde{G}}^{\delta }(\mathbf{x})-{\tilde{G}}^{\ast }(\mathbf{x});$(33) ) and (34(34) $\Vert \mathcal{E}{\Vert }_{H^3(\mathrm{\Omega }{)}^N}<\delta \mathrm{as}\delta \to 0^{+}$(34) ). Then, we have the estimate (35) $\Vert U^{\delta }-U^{\ast }{\Vert }_{H^1(\mathrm{\Omega }{)}^N}^2\le C({\delta }^2+ϵ\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2),$(35) where C is a constant that depends only on Ω, $\Vert A{\Vert }_{C^1(\overline{\mathrm{\Omega }})}$ and μ.

Proof.

Since $U^{\delta }$ is the minimizer of $J_{ϵ}$, by the variational principle, we have (36) $⟨\mathcal{L}{U}^{\delta }-S{U}^{\delta },\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{⟩}_{L^2(\mathrm{\Omega }{)}^N}+ϵ⟨U^{\delta },\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}=0$(36) for all test functions Φ in the space $H_0^3(\mathrm{\Omega }{)}^N=\{\varphi \in H^3(\mathrm{\Omega }{)}^N:\varphi =A\mathrm{\nabla }\varphi \cdot \nu =0\mathrm{on}\mathrm{\partial }\mathrm{\Omega }\}.$ Since $\mathcal{L}{U}^{\ast }-S{U}^{\ast }=0,$ we can deduce from (36(36) $⟨\mathcal{L}{U}^{\delta }-S{U}^{\delta },\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{⟩}_{L^2(\mathrm{\Omega }{)}^N}+ϵ⟨U^{\delta },\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}=0$(36) ) that $\begin{array}{rl}& ⟨\mathcal{L}(U^{\delta }-U^{\ast })-S(U^{\delta }-U^{\ast }),\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{⟩}_{L^2(\mathrm{\Omega }{)}^N}+ϵ⟨U^{\delta }-U^{\ast },\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}=\\ & -ϵ⟨U^{\ast },\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}.\end{array}$ Plugging the test function (37) $\mathrm{\Phi }=U^{\delta }-U^{\ast }-\mathcal{E}\in H_0^3(\mathrm{\Omega })$(37) into the identity above, we have $\begin{array}{rl}\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2+ϵ\Vert \mathrm{\Phi }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2& =-⟨\mathcal{L}\mathcal{E}-S\mathcal{E},\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{⟩}_{L^2(\mathrm{\Omega }{)}^N}^2\\ & -ϵ⟨\mathcal{E},\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}-ϵ⟨U^{\ast },\mathrm{\Phi }{⟩}_{H^3(\mathrm{\Omega }{)}^N}.\end{array}$ Applying the Cauchy–Schwartz inequality and removing lower order terms, we obtain (38) $\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2+ϵ\Vert \mathrm{\Phi }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2\le C({\delta }^2+ϵ\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2).$(38) Recall from (3(3) $\mathcal{L}v=\mathrm{Div}(A\mathrm{\nabla }v)+\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }v(\mathbf{x})+c(\mathbf{x})v$(3) ) that $\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2=\Vert \mathrm{Div}(A(\mathbf{x})\mathrm{\nabla }\mathrm{\Phi }+\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }\mathrm{\Phi }+(c(\mathbf{x})\mathrm{Id}-S)\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2.$ Recall the function ψ in (24(24) $\psi (\mathbf{x})=x_1+\frac{1}{2X^2}\sum _{i=1}^d{x}_i^2+\alpha ,0<\alpha <1/2.$(24) ). Fix $\lambda >{\sigma }_0$ and $p>{\sigma }_1$ as in Lemma 3.1. Set $M=max\{e^{2\lambda {\psi }^{-p}(\mathbf{x})}:\mathbf{x}\in \overline{\mathrm{\Omega }}\}\mathrm{and}m=min\{e^{2\lambda {\psi }^{-p}(\mathbf{x})}:\mathbf{x}\in \overline{\mathrm{\Omega }}\}.$ We have $\begin{array}{rl}M\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2& \ge \frac{1}{2}\Vert e^{\lambda {\psi }^{-p}(\mathbf{x})}\mathrm{Div}(A(\mathbf{x})\mathrm{\nabla }\mathrm{\Phi })\\ & +e^{\lambda {\psi }^{-p}(\mathbf{x})}(\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }\mathrm{\Phi }+(c(\mathbf{x})\mathrm{Id}-S)\mathrm{\Phi }){\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2.\end{array}$ Using the inequality $(x-y{)}^2\ge \frac{1}{2}{x}^2-y^2$, we have $\begin{array}{rl}M\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2& \ge \frac{1}{2}\Vert e^{\lambda {\psi }^{-p}(\mathbf{x})}\mathrm{Div}(A(\mathbf{x})\mathrm{\nabla }\mathrm{\Phi }){\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2\\ & -\Vert e^{\lambda {\psi }^{-p}(\mathbf{x})}(\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }\mathrm{\Phi }+(c(\mathbf{x})\mathrm{Id}-S)\mathrm{\Phi }){\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2.\end{array}$ Hence, thus, by (29(29) $\begin{array}{rl}& \lambda p{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }u{|}^2\mathrm{d}\mathbf{x}+{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|u{|}^2\mathrm{d}\mathbf{x}\\ & \le C{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{Div}(A\mathrm{\nabla }u){|}^2\mathrm{d}\mathbf{x}.\end{array}$(29) ), (39) $\begin{array}{rl}M\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2& \ge C\lambda p{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }\mathrm{\Phi }{|}^2\mathrm{d}\mathbf{x}+C{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\Phi }{|}^2\mathrm{d}\mathbf{x}\\ & -\Vert e^{\lambda {\psi }^{-p}(\mathbf{x})}(\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }\mathrm{\Phi }+(c(\mathbf{x})\mathrm{Id}-S)\mathrm{\Phi }){\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2.\end{array}$(39) Now, fixing $\lambda >{\sigma }_0$ large, we obtain from (39(39) $\begin{array}{rl}M\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2& \ge C\lambda p{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\nabla }\mathrm{\Phi }{|}^2\mathrm{d}\mathbf{x}+C{\lambda }^3{p}^4{\int }_{\mathrm{\Omega }}{e}^{2\lambda {\psi }^{-p}(\mathbf{x})}|\mathrm{\Phi }{|}^2\mathrm{d}\mathbf{x}\\ & -\Vert e^{\lambda {\psi }^{-p}(\mathbf{x})}(\mathbf{b}(\mathbf{x})\cdot \mathrm{\nabla }\mathrm{\Phi }+(c(\mathbf{x})\mathrm{Id}-S)\mathrm{\Phi }){\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2.\end{array}$(39) ) that (40) $M|\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2\ge Cm\Vert \mathrm{\Phi }{\Vert }_{H^1(\mathrm{\Omega }{)}^N}.$(40) Here, we have used the boundedness of $\mathbf{b}$ and c in Ω. Combining (37(37) $\mathrm{\Phi }=U^{\delta }-U^{\ast }-\mathcal{E}\in H_0^3(\mathrm{\Omega })$(37) ), (38(38) $\Vert \mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2+ϵ\Vert \mathrm{\Phi }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2\le C({\delta }^2+ϵ\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2).$(38) ) and (40(40) $M|\mathcal{L}\mathrm{\Phi }-S\mathrm{\Phi }{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2\ge Cm\Vert \mathrm{\Phi }{\Vert }_{H^1(\mathrm{\Omega }{)}^N}.$(40) ) gives $\Vert U^{\delta }-U^{\ast }-\mathcal{E}{\Vert }_{H^1(\mathrm{\Omega }{)}^N}^2\le C({\delta }^2+ϵ\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2).$ This and the assumption $\Vert \mathcal{E}{\Vert }_{H^3(\mathrm{\Omega }{)}^N}\le C\delta$ imply inequality (35(35) $\Vert U^{\delta }-U^{\ast }{\Vert }_{H^1(\mathrm{\Omega }{)}^N}^2\le C({\delta }^2+ϵ\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}^2),$(35) ).

Corollary 4.2

Let $f^{\ast }(\mathbf{x})=u^{\ast }(\mathbf{x},0)$ and $f^{\delta }(\mathbf{x})=u^{\delta }(\mathbf{x},0)$ where $u^{\ast }(\mathbf{x},t)$ and $u^{\delta }(\mathbf{x},t)$ are computed from $U^{\ast }(\mathbf{x})$ and $U^{\delta }(\mathbf{x})$ via (8(8) $u(\mathbf{x},t)=\sum _{n=1}^{\mathrm{\infty }}{u}_n(\mathbf{x}){\mathrm{\Psi }}_n(t),$(8) ) and (17(17) $U(\mathbf{x})=(u_1(\mathbf{x}),u_2(\mathbf{x}),\dots ,u_N(\mathbf{x}){)}^T.$(17) ). Then, by the trace theory $\Vert f^{\delta }-f^{\ast }{\Vert }_{L^2(\mathrm{\Omega })}\le C(\delta +\sqrt{ϵ}\Vert U^{\ast }{\Vert }_{H^3(\mathrm{\Omega }{)}^N}).$

5. Numerical illustrations

We numerically test our method when d=2. The domain Ω is the square $(-R,R{)}^2$. In this section, we write $\mathbf{x}=(x,y)$. For the coefficients of the governing equation, we choose, for simplicity, $A(\mathbf{x})=\mathrm{Id}$ and $\mathbf{b}(\mathbf{x})=0.$ The function c is set as $c(x,y)=(3(1-x{)}^2{e}^{-x^2-(y+1{)}^2}-10(x/5-x^3-y^5)e^{-x^2-y^2}-1/3e^{-(x+1{)}^2-y^2})/10$ which is a scale of the ‘peaks’ function in Matlab. The graph of c is displayed in Figure 2.

Figure 2. The true function $c(\mathbf{x})$ used for all numerical examples in this section.

Define a grid of points in Ω $\mathcal{G}=\{(x_i,y_j)=(-R+(i-1)d_{\mathbf{x}},-R+(j-1)d\mathbf{x}):1\le i,j\le N_{\mathbf{x}}+1\},$ where $N_{\mathbf{x}}=80$ and $d_{\mathbf{x}}=2R/N_{\mathbf{x}}.$ For the time variable, we choose T=4. Define a uniform partition of $[0,T]$ as $0=t_1<t_2<\cdots <t_{N_T+1}=T$ with step size $d_t=T/N_T$. In our tests, $N_T=250.$ The forward problem is solved by finite difference method in the implicit scheme. Denote by $u^{\ast }$ the solution of the forward problem. The data are given by $F(\mathbf{x},t)=u(\mathbf{x},t)(1+\delta (2\mathrm{rand}-1)),G(\mathbf{x},t)={\mathrm{\partial }}_{\nu }u(\mathbf{x},t)(1+\delta (2\mathrm{rand}-1))$ for $(\mathbf{x},t)\in \mathrm{\partial }\mathrm{\Omega }\times [0,T]$ where $\mathrm{rand}$ is the uniformly distributed random function taking value in $[0,1]$ and δ is the noise level. The noise level δ is given in each numerical tests.

5.1. The implementation for Algorithm 1

The main part of this section is to compute the minimizer U of $J_{ϵ}$ subject to the constraints (19(19) $\begin{array}{rl}U(\mathbf{x})& =\tilde{F}(\mathbf{x})={({\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}F(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(19) ) and (20(20) $\begin{array}{rl}A\mathrm{\nabla }U(\mathbf{x})\cdot \nu & =\tilde{G}(\mathbf{x})={({\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_1(t)\mathrm{d}t,\dots ,{\int }_0^{T}G(\mathbf{x},t){\mathrm{\Psi }}_N(t)\mathrm{d}t)}^T\end{array}$(20) ). The ‘cut-off’ number N is set to be 30, see Remark 2.2 for this choice of N. To construct the orthonormal basis $\{{\mathrm{\Psi }}_n{\}}_{n=1}^N,$ for each $n\in \{1,\dots ,N\}$, we identify the function ${\mathrm{\Phi }}_n$, defined in (37(37) $\mathrm{\Phi }=U^{\delta }-U^{\ast }-\mathcal{E}\in H_0^3(\mathrm{\Omega })$(37) ), by the $N_T+1$ dimensional vector $({\mathrm{\Phi }}_n(t_1),\dots ,{\mathrm{\Phi }}_n(t_{N_T+1}))$. Then apply the Gram–Schmidt orthogonalization process for the set $\{({\mathrm{\Phi }}_n(t_1),\dots ,{\mathrm{\Phi }}_n(t_{N_T+1})){\}}_{n=1}^N$ in the $N_T+1$ dimensional Euclidian space. In other words, we construct $\{{\mathrm{\Psi }}_n{\}}_{n=1}^N$ in the finite difference scheme. The discretized version of $U(\mathbf{x})=(u_1(\mathbf{x}),\dots ,u_N(\mathbf{x}){)}^T$, $\mathbf{x}\in \mathrm{\Omega }$ is $(u_1(x_i,y_j),\dots ,u_N(x_i,y_j){)}_{i,j=1}^{N_{\mathbf{x}}+1}.$ Hence, $J_{ϵ}(U)$, see (21(21) $J_{ϵ}(U)={\int }_{\mathrm{\Omega }}|\mathcal{L}U(\mathbf{x})-SU(\mathbf{x}){|}^2\mathrm{d}\mathbf{x}+ϵ\Vert U{\Vert }_{H^3(\mathrm{\Omega })}^2$(21) ), is approximated by (41) $\begin{array}{rl}{J}_{ϵ}(U)& =d_{\mathbf{x}}^2\sum _{i,i=2}^{N_{\mathbf{x}}}\sum _{m=1}^N|\frac{u_m(x_{i+1},y_j)+u_m(x_{i-1},y_j)+u_m(x_i,y_{j+1})+u_m(x_i,y_{j-1})-4u_m(x_i,y_j)}{d_{\mathbf{x}}^2}\\ & {+c(x_i,y_j)u_m(x_i,y_j)-\sum _{n=1}^N{s}_{mn}{u}_n(x_i,y_j)|}^2+ϵd_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|u_m(x_i,y_j){|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_{i+1},y_j)-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_i,y_{j+1})-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2.\end{array}$(41) Here, we slightly change the $H^3$ norm of the regularity term to the $H^1$ norm. This makes the computational codes less heavy. The numerical results with this change are still acceptable. We also modify the regularized parameter of the term $\Vert \mathrm{\nabla }U{\Vert }_{L^2(\mathrm{\Omega }{)}^N}$ to be ${ϵ}^2$, instead of ε, since we observe that the obtained numerical results are more accurate with this modification. To numerically prove this, we solve the inverse problem when the function $f_{\mathrm{true}}$ is given in Test 1 in Section 5.2 in two cases: with and without this modification and then compare the corresponding outputs. The results are displayed in Figure 3. It is clear from Figure 3 that the modification above provides better numerical results.

Figure 3. Test 1. The comparison of the reconstruction of the function f with and without the modification for the regularized parameter. It is evident that the numerical result in (b) is significantly better than that in (c) in both reconstructed shape and computed value. (a) The function $f_{\mathrm{true}}$. (b) $f_{\mathrm{comp}}$ computed using the regularization term ${ϵ}^2\Vert \mathrm{\nabla }U{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2$ when $\delta =50\mathrm{\%}$. (c) $f_{\mathrm{comp}}$ computed using the regularization term $ϵ\Vert \mathrm{\nabla }U{\Vert }_{L^2(\mathrm{\Omega }{)}^N}^2$ when $\delta =50\mathrm{\%}$.

The expression in (41(41) $\begin{array}{rl}{J}_{ϵ}(U)& =d_{\mathbf{x}}^2\sum _{i,i=2}^{N_{\mathbf{x}}}\sum _{m=1}^N|\frac{u_m(x_{i+1},y_j)+u_m(x_{i-1},y_j)+u_m(x_i,y_{j+1})+u_m(x_i,y_{j-1})-4u_m(x_i,y_j)}{d_{\mathbf{x}}^2}\\ & {+c(x_i,y_j)u_m(x_i,y_j)-\sum _{n=1}^N{s}_{mn}{u}_n(x_i,y_j)|}^2+ϵd_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|u_m(x_i,y_j){|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_{i+1},y_j)-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_i,y_{j+1})-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2.\end{array}$(41) ) is simplified as follows: (42) $\begin{array}{rl}{J}_{ϵ}(U)& =d_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|\sum _{n=1}^N[{\delta }_{mn}(\frac{-4}{d_{\mathbf{x}}^2}+c(x_i,y_j))-s_{mn}]u_n(x_i,y_j)+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_{i+1},y_j)\\ & {+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_{i-1},y_j)+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_i,y_{j+1})+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_i,y_{j-1})|}^2\\ & +ϵd_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|u_m(x_i,y_j){|}^2+{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_{i+1},y_j)-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_i,y_{j+1})-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2.\end{array}$(42) Here, we use the Kronecker number ${\delta }_{mn}$ for the convenience of writing the computational codes. We next identify $\{u_n(x_i,y_j):1\le i,j\le N_{\mathbf{x}}+1,1\le n\le N\}$ with the $(N_{\mathbf{x}}+1{)}^{2}N$ dimensional vector $\mathfrak{U}=({\mathfrak{u}}_{\mathfrak{i}}{)}_{\mathfrak{i}=1}^{(N_{\mathbf{x}}+1{)}^{2}N}$ according to the rule ${\mathfrak{u}}_{\mathfrak{i}}=u_n(x_i,y_j)$ where the index $\mathfrak{i}$ is $\mathfrak{i}=(i-1)(N_{\mathbf{x}}+1)N+(j-1)N+n,1\le i,j\le N_{\mathbf{x}}+1,1\le n\le N.$ Then, with this notation, $J_{ϵ}(U)$ in (42(42) $\begin{array}{rl}{J}_{ϵ}(U)& =d_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|\sum _{n=1}^N[{\delta }_{mn}(\frac{-4}{d_{\mathbf{x}}^2}+c(x_i,y_j))-s_{mn}]u_n(x_i,y_j)+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_{i+1},y_j)\\ & {+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_{i-1},y_j)+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_i,y_{j+1})+\frac{{\delta }_{mn}}{d_{\mathbf{x}}^2}{u}_n(x_i,y_{j-1})|}^2\\ & +ϵd_{\mathbf{x}}^2\sum _{i,j=2}^{N_x}\sum _{m=1}^N|u_m(x_i,y_j){|}^2+{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_{i+1},y_j)-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2\\ & +{ϵ}^2{d}_{\mathbf{x}}^2\sum _{i,j=2}^{N_{\mathbf{x}}}\sum _{m=1}^N{\left|\frac{u_m(x_i,y_{j+1})-u_m(x_i,y_j)}{d_{\mathbf{x}}}\right|}^2.\end{array}$(42) ) is rewritten as ${\mathfrak{J}}_{ϵ}(\mathfrak{U})=d_{\mathbf{x}}^2|\mathfrak{L}\mathfrak{U}{|}^2++ϵd_{\mathbf{x}}^2|\mathfrak{U}{|}^2+ϵd_{\mathbf{x}}^2|D_x\mathfrak{U}{|}^2+ϵd_{\mathbf{x}}^2|D_y\mathfrak{U}{|}^2.$ The $(N_{\mathbf{x}}+1{)}^{2}N\times (N_{\mathbf{x}}+1{)}^{2}N$ matrices $\mathfrak{L}$, $D_x$ and $D_y$ are as follows.

1. Define the matrix $\mathfrak{L}.$ For $\mathfrak{i}=(i-1)(N_{\mathbf{x}}+1)N+(j-1)N+m$, for some $2\le i,j\le N_x$, the $\mathfrak{i}\mathfrak{j}\mathrm{th}$ entry of $\mathfrak{L}$ is

   1. ${\delta }_{mn}(-4/d_{\mathbf{x}}^2+c(x_i,y_j))-s_{mn}$ if $\mathfrak{j}=(i-1)(N_{\mathbf{x}}+1)N+(j-1)N+n,$

   2. ${\delta }_{mn}/d_x^2$ if $\mathfrak{j}=(i\pm 1-1)(N_{\mathbf{x}}+1)N+(j-1)N+n$ or $\mathfrak{j}=(i-1)(N_{\mathbf{x}}+1)N+(j\pm 1-1)N+n$,

   3. 0 otherwise.

2. Define the matrix $D_x$. For $\mathfrak{i}=(i-1)(N_{\mathbf{x}}+1)N+(j-1)N+m$, for some $2\le i,j\le N_x$, the $\mathfrak{i}\mathfrak{j}\mathrm{th}$ entry of $D_x$ is

   1. $1/d_{\mathbf{x}}$ if $\mathfrak{j}=(i+1-1)(N_{\mathbf{x}}+1)N+(j-1)N+m$,

   2. $-1/d_{\mathbf{x}}$ if $\mathfrak{j}=\mathfrak{i}$,

   3. 0 otherwise.

3. Define the matrix $D_x$. For $\mathfrak{i}=(i-1)(N_{\mathbf{x}}+1)N+(j-1)N+m$, for some $2\le i,j\le N_x$, the $\mathfrak{i}\mathfrak{j}\mathrm{th}$ entry of $D_x$ is

   1. $1/d_{\mathbf{x}}$ if $\mathfrak{j}=(i-1)(N_{\mathbf{x}}+1)N+(j+1-1)N+m$,

   2. $-1/d_{\mathbf{x}}$ if $\mathfrak{j}=\mathfrak{i}$,

   3. 0 otherwise.

Remark 5.1

The values of the parameters

As mentioned, we take N=30, $N_{\mathbf{x}}=80$, $N_T=250$, R=2. The regularized parameter $ϵ={10}^{-7}$. These values of parameters are used for all tests in Section 5.2.

5.2. Tests

We perform four (4) numerical examples in this paper. These examples with high levels of noise show the strength of our method. We will also compare the reconstructed maximum values of the reconstructed functions and the true ones. Below, $f_{\mathrm{true}}$ and $f_{\mathrm{comp}}$ are, respectively, the true source function and the reconstructed one due to Algorithm 1 with the parameters in Section 5.1.

1. Test 1. The case of one inclusion. The function $f_{\mathrm{true}}$ is a smooth function supported in a disk with radius 1 centred at the origin. More precisely, $f_{\mathrm{true}}(\mathbf{x})=\{\begin{array}{ll}{\displaystyle e^{-(1/(1-|\mathbf{x}{|}^2))+1}}& \mathrm{if}|\mathbf{x}|<1,\\ 0& \mathrm{otherwise}.\end{array}$ Figure 4 displays the functions $f_{\mathrm{true}}$ and $f_{\mathrm{comp}}$. Table 1 show the reconstructed value of the function $f_{\mathrm{comp}}$ and the relative error. The noise levels are $\delta =0\mathrm{\%}$, $25\mathrm{\%}$, $50\mathrm{\%}$, $75\mathrm{\%}$ and $100\mathrm{\%}.$

   Table 1. Test 1. Correct and computed maximal values of source functions.

   Noise level (%)	$max{f}_{\mathrm{true}}$	$max{f}_{\mathrm{comp}}$	${\mathrm{error}}_{\mathrm{rel}}$ (%)	${\mathrm{pos}}_{\mathrm{true}}$	${\mathrm{pos}}_{\mathrm{comp}}$	$\mathrm{di}{\mathrm{s}}_{\mathrm{err}}$
   0	1	0.99	1.0	$(0,0)$	$(0,0)$	0
   25	1	0.96	4.0	$(0,0)$	$(-0.05,0)$	0.05
   50	1	1.21	21.0	$(0,0)$	$(-0.05,0.1)$	0.11
   75	1	1.01	1.0	$(0,0)$	$(0.2,-0.1)$	0.22
   100	1	1.53	53.0	$(0,0)$	$(0.25,0.1)$	0.27

   Notes: ${\mathrm{error}}_{\mathrm{rel}}$ denotes the relative error of the reconstructed maximal value. ${\mathrm{pos}}_{\mathrm{true}}$ is the true position of the inclusion; i.e. the maximizer of $f_{\mathrm{true}}$. ${\mathrm{pos}}_{\mathrm{comp}}$ is the computed position of the inclusion. ${\mathrm{dis}}_{\mathrm{err}}$ is the absolute error of the reconstructed positions.

   Figure 4. Test 1. The true and computed source functions. Our method still works well when $\delta =100\mathrm{\%}.$ It is shown in (e) that the reconstructed value of $f_{\mathrm{comp}}$ with $\delta =75\mathrm{\%}$ is quite accurate, even better than in (d), but in contrast, the reconstructed shape starts to break down. (a) The function $f_{\mathrm{true}}$. (b) $f_{\mathrm{comp}}$, $\delta =0\mathrm{\%}$. (c) $f_{\mathrm{comp}}$, $\delta =25\mathrm{\%}$. (d) $f_{\mathrm{comp}}$, $\delta =50\mathrm{\%}$. (e) $f_{\mathrm{comp}}$, $\delta =75\mathrm{\%}$. (f) $f_{\mathrm{comp}}$, $\delta =100\mathrm{\%}$.

   

   It is evident that our method is robust for Test 1 in the sense that the reconstructed maximal value of the function f and the reconstructed shape and position of the inclusion are quite accurate.

2. Test 2. The case of two inclusions. The function $f_{\mathrm{true}}$ is a smooth function supported in two disks with radius r=0.8 centred at ${\mathbf{x}}_1=(-1,0)$ and ${\mathbf{x}}_2=(1,0)$ respectively. The function $f_{\mathrm{true}}$ is given by the formula $f(\mathbf{x})=\{\begin{array}{ll}{\displaystyle e^{-(r^2/(r^2-|\mathbf{x}-{\mathbf{x}}_1{|}^2))+1}}& \mathrm{if}|\mathbf{x}-{\mathbf{x}}_1|<r,\\ {\displaystyle e^{-(r^2/(r^2-|\mathbf{x}-{\mathbf{x}}_2{|}^2))+1}}& \mathrm{if}|\mathbf{x}-{\mathbf{x}}_2|<r,\\ 0& \mathrm{otherwise}.\end{array}$ Figure 5 displays the functions $f_{\mathrm{true}}$ and $f_{\mathrm{comp}}$. Table 2 show the reconstructed value of the function $f_{\mathrm{comp}}$ and the relative error. The noise levels are $\delta =0\mathrm{\%}$, $25\mathrm{\%}$, $50\mathrm{\%}$, $75\mathrm{\%}$ and $100\mathrm{\%}.$ The reconstruction in Test 2 is good. In this test, the reconstruct breaks down when the noise level is $75\mathrm{\%}$ although we are able to detect the inclusions with higher noise levels.

3. Test 3. The case of non-inclusion and nonsmooth function. The function $f_{\mathrm{true}}$ is the characteristic function of the letter Y. Figure 6 displays the functions $f_{\mathrm{true}}$ and $f_{\mathrm{comp}}$. The noise levels are $\delta =10\mathrm{\%}$ and $15\mathrm{\%}.$

   Table 2. Test 2. Correct and computed maximal values of the inclusions.

   Noise level (%)	Inclusion	$\underset{\mathrm{inc},\mathrm{true}}{max}$	$\underset{\mathrm{inc},\mathrm{comp}}{max}$	${\mathrm{error}}_{\mathrm{rel}}$ (%)	${\mathrm{pos}}_{\mathrm{true}}$	${\mathrm{pos}}_{\mathrm{comp}}$	$\mathrm{di}{\mathrm{s}}_{\mathrm{err}}$
   0	Left	1	0.96	4.0	$(-1,0)$	$(-0.85,0)$	0.15
   0	Right	1	0.98	2.0	$(1,0)$	$(0.85,0)$	0.15
   25	Left	1	1.11	11	$(-1,0)$	$(-0.85,0)$	0.15
   25	Right	1	0.96	4	$(1,0)$	$(0.9,0.1)$	0.14
   50	Left	1	0.61	39	$(-1,0)$	$(-0.9,0.25)$	0.27
   50	Right	1	1.01	1	$(1,0)$	$(0.85,-0.2)$	0.25
   75	Left	1	0.84	26	$(-1,0)$	$(-1,0.1)$	0.1
   75	Right	1	1.82	82	$(1,0)$	$(0.8,0)$	0.2
   100	Left	1	1.1	10	$(-1,0)$	$(-0.9,0.3)$	0.32
   100	Right	1	1.58	58	$(1,0)$	$(0.05,0.8)$	0.21

   Notes: $\underset{\mathrm{inc},\mathrm{true}}{max}$ and $\underset{\mathrm{inc},\mathrm{comp}}{max}$ are the true and computed, respectively, maximal values of the source in an inclusion. ${\mathrm{error}}_{\mathrm{rel}}$ denotes the relative error of the reconstructed maximal value. ${\mathrm{pos}}_{\mathrm{true}}$ is the true position of the inclusion; i.e. the maximizer of $f_{\mathrm{true}}$. ${\mathrm{pos}}_{\mathrm{comp}}$ is the computed position of the inclusion. ${\mathrm{dis}}_{\mathrm{err}}$ is the absolute error of the reconstructed positions.

   Figure 5. Test 2. The true and computed source functions. The reconstruction of the two inclusions are not symmetric probably because the true function c, see Figure 2 for its graph, is negative on the left and positive on the right. However, both inclusions can be seen when the noise level goes up to $100\mathrm{\%}$. (a) The function $f_{\mathrm{true}}$. (b) $f_{\mathrm{comp}}$, $\delta =0\mathrm{\%}$. (c) $f_{\mathrm{comp}}$, $\delta =25\mathrm{\%}$. (d) $f_{\mathrm{comp}}$, $\delta =50\mathrm{\%}$. (e) $f_{\mathrm{comp}}$, $\delta =75\mathrm{\%}$. (f) $f_{\mathrm{comp}}$, $\delta =100\mathrm{\%}$.

   

   Figure 6. Test 3. The true and computed source functions. The letter Y can be detected well in this case. The true maximal value of $f_{\mathrm{true}}$ is 1. The computed maximal value of $f_{\mathrm{comp}}$ when $\delta =10\mathrm{\%}$ is 1.09 (relative error 9%). The computed maximal value of $f_{\mathrm{comp}}$ when $\delta =15\mathrm{\%}$ is 1.15 (relative error 15%). (a) The function $f_{\mathrm{true}}$. (b) $f_{\mathrm{comp}}$, $\delta =10\mathrm{\%}$. (c) $f_{\mathrm{comp}}$, $\delta =15\mathrm{\%}$.

   

   We can reconstruct the letter Y and the reconstructed maximal of $f_{\mathrm{comp}}$ is good when $\delta =10\mathrm{\%}$ but the error is large when the noise level reaches $15\mathrm{\%}.$

4. Test 4. The case of non-inclusion and nonsmooth function. The function $f_{\mathrm{true}}$ is the characteristic function of the letter λ. Figure 7 displays the functions $f_{\mathrm{true}}$ and $f_{\mathrm{comp}}$. The noise levels are $\delta =10\mathrm{\%}$ and $15\mathrm{\%}.$

   Figure 7. Test 4. The true and computed source functions. The reconstruction of λ is acceptable. The true maximal value of $f_{\mathrm{true}}$ is 1. The computed maximal value of $f_{\mathrm{comp}}$ when $\delta =10\mathrm{\%}$ is 1.16 (relative error 16%). The computed maximal value of $f_{\mathrm{comp}}$ when $\delta =15\mathrm{\%}$ is 1.11 (relative error 11%). (a) The function $f_{\mathrm{true}}$. (b) $f_{\mathrm{comp}}$, $\delta =10\mathrm{\%}$. (c) $f_{\mathrm{comp}}$, $\delta =15\mathrm{\%}$.

   

   The image of λ in Test 4 is acceptable. The reconstructed maximal value in Figure 7(c) is better than that in Figure 7(b) but the reconstruction of λ in Figure 7(c) is not as good as that in Figure 7(b).

6. Concluding remarks

In this paper, we have solved the problem of reconstructing the initial condition of solution to a general class of parabolic equation from the measurement of lateral Cauchy data. The main points of the method is derive an approximate model by a truncation of the Fourier series with respect to a special basis. We solved the approximation model by the quasi-reversibility method. The convergence of this method when the noise tends to 0 was proved. More importantly, numerical examples show that our method is robust when proving accurate reconstructions of the unknown source function from highly noisy data.

Although our method leads to good numerical results, it has a drawback. The proof of the ‘convergence’ of the system (18(18) $\mathcal{L}U(\mathbf{x})=SU(\mathbf{x})\mathrm{for}\mathbf{x}\in \mathrm{\Omega },$(18) ) as $N\to \mathrm{\infty }$ is challenging and is omitted in this paper. We refer the reader to [29, Section 4] for an alternative approach to solve Problem 1.1 by which we can avoid this non-rigorousness. This method is based on the Carleman estimate for parabolic operators. However, in this case we can determine a ‘near’ initial condition for the function $u(\mathbf{x},t)$. That means, we can recover the function $u(\mathbf{x},ϵ)$ where ε is any small number. Implementation for the method in [29, Section 4] is valuable. We reserve it for a future research.

Acknowledgments

The authors are grateful to Michael V. Klibanov for many fruitful discussions.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Loc Hoang Nguyen  http://orcid.org/0000-0002-0172-8816

Additional information

Funding

The work of Nguyen was supported by the US Army Research Laboratory and US Army Research Office grant W911NF-19-1-0044. In addition, the effort of Nguyen and Li was supported by research funds FRG 111172 provided by The University of North Carolina at Charlotte.