A new algorithm to determine the creation or depletion term of parabolic equations from boundary measurements

TL;DR

This paper introduces a novel numerical method for determining creation or depletion coefficients in parabolic equations from boundary data, avoiding the need for prior knowledge of the true solution and demonstrating effectiveness through numerical tests.

Contribution

A new approach that solves nonlinear inverse problems for parabolic equations without requiring an initial guess of the true coefficient.

Findings

Method successfully recovers coefficients from boundary measurements

Numerical experiments confirm robustness and accuracy

Avoids reliance on advanced prior knowledge of solutions

Abstract

We propose a robust numerical method to find the coefficient of the creation or depletion term of parabolic equations from the measurement of the lateral Cauchy information of their solutions. Most papers in the field study this nonlinear and severely ill-posed problem using optimal control. The main drawback of this widely used approach is the need of some advanced knowledge of the true solution. In this paper, we propose a new method that opens a door to solve nonlinear inverse problems for parabolic equations without any initial guess of the true coefficient. This claim is confirmed numerically. The key point of the method is to derive a system of nonlinear elliptic equations for the Fourier coefficients of the solution to the governing equation with respect to a special basis of L^2. We then solve this system by a predictor-corrector process, in which our computation to obtain the first and second predictors is effective. The desired solution to the inverse problem under consideration follows.

Equations

$$A(\mathbf{x})\xi \cdot \xi \ge \mu |\xi {|}^2\text{mbox}forall\mathbf{x}\in {\mathbb{R}}^d,\xi \in {\mathbb{R}}^d,$$

$$Lw(\mathbf{x})=\mathrm{div}(A(\mathbf{x})\nabla w(\mathbf{x}))+\mathbf{b}(\mathbf{x})\cdot \nabla w(\mathbf{x})\text{mbox}forallw\in H^2({\mathbb{R}}^d),\mathbf{x}\in {\mathbb{R}}^d.$$

$$\left\{\begin{array}[]{rcll}u_{t}({\bf x},t)&=&Lu({\bf x},t)+c({\bf x})u({\bf x},t)&{\bf x}\in\mathbb{R}^{d},t\in(0,\infty),\\ u({\bf x},0)&=&f({\bf x})&{\bf x}\in\mathbb{R}^{d}.\end{array}\right.$$

$$F(\mathbf{x},t)=u(\mathbf{x},t)\text{mbox}andG(\mathbf{x},t)=A(\mathbf{x})\nabla u(\mathbf{x},t)\cdot \nu (\mathbf{x})$$

$$v(\mathbf{x},t)=u_t(\mathbf{x},t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$v_t(\mathbf{x},t)=Lv(\mathbf{x},t)+c(\mathbf{x})v(\mathbf{x},t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$v(\mathbf{x},0)=u_t(\mathbf{x},0)=Lf(\mathbf{x})+c(\mathbf{x})f(\mathbf{x}).$$

$$c(\mathbf{x})=\frac{v(\mathbf{x},0)-Lf(\mathbf{x})}{f(\mathbf{x})}\text{mbox}forall\mathbf{x}\in \Omega .$$

$$v_t(\mathbf{x},t)=Lv(\mathbf{x},t)-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}v(\mathbf{x},t)+\frac{v(\mathbf{x},0)}{f(\mathbf{x})}v(\mathbf{x},t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$s_{mn}=\int_{0}^{T}\Psi_{n}^{\prime}(t)\Psi_{m}(t)dt=\left\{\begin{array}[]{ll}1&\mbox{if }m=n,\\ 0&\mbox{if }n<m.\end{array}\right.$$

$$v_n(\mathbf{x})={\int }_0^{T}v(\mathbf{x},t){\Psi }_n(t)dt\text{mbox}forall\mathbf{x}\in \Omega .$$

$$v(\mathbf{x},t)=\sum _{n=1}^{\infty }{v}_n(\mathbf{x}){\Psi }_n(t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$v(\mathbf{x},t)=\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$v(\mathbf{x},0)=\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(0)\text{mbox}forall\mathbf{x}\in \Omega$$

$$v_t(\mathbf{x},t)=\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n^{\prime }(t)\text{mbox}forall(\mathbf{x},t)\in {\Omega }_T.$$

$$\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n^{\prime }(t)=\sum _{n=1}^{N}L{v}_n(\mathbf{x}){\Psi }_n(t)-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(t)+\sum _{n,l=1}^N\frac{v_n(\mathbf{x}){\Psi }_n(0)}{f(\mathbf{x})}{v}_l(\mathbf{x}){\Psi }_l(t)$$

$${\int }_0^T{\Psi }_n(t){\Psi }_m(t)dt={\delta }_{m-n},$$

$$L{v}_m(\mathbf{x})-\sum _{n=1}^N{s}_{mn}{v}_n(\mathbf{x})-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}{v}_m(\mathbf{x})+\sum _{n=1}^N\frac{{\Psi }_n(0)}{f(\mathbf{x})}{v}_n(\mathbf{x})v_m(\mathbf{x})=0$$

$$\begin{array}[]{c}\displaystyle v_{m}({\bf x})=F_{m}({\bf x})=\int_{0}^{T}F_{t}({\bf x},t)\Psi_{m}(t)dt,\\ \displaystyle\partial_{\nu}v_{m}({\bf x})=G_{m}({\bf x})=\int_{0}^{T}G_{t}({\bf x},t)\Psi_{m}(t)dt,\end{array}$$

$$F^{\delta }=F^{\ast }+\delta \mathrm{rand}{G}^{\delta }=G^{\ast }+\delta \mathrm{rand}.$$

$${\int }_{\Omega }|\mathrm{div}(A\nabla u){|}^2\exp \left[2\lambda {\left(x_d+b\right)}^p\right]d\mathbf{x}\ge C\lambda {\int }_{\Omega }\left[|\nabla u{|}^2+{\lambda }^2{u}^2\right]\exp \left[2\lambda {\left(x_d+b\right)}^p\right]d\mathbf{x},$$

$$L{u}_m(\mathbf{x})-\sum _{n=1}^N{s}_{mn}{u}_n(\mathbf{x})-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}{u}_m(\mathbf{x})+\sum _{n=1}^N\frac{{\Psi }_n(0)}{f(\mathbf{x})}{u}_n(\mathbf{x})u_m(\mathbf{x})=0$$

$$L{v}_m(\mathbf{x})-\sum _{n=1}^N{s}_{mn}{v}_n(\mathbf{x})-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}{v}_m(\mathbf{x})+\sum _{n=1}^N\frac{{\Psi }_n(0)}{f(\mathbf{x})}{v}_n(\mathbf{x})v_m(\mathbf{x})=0$$

$$\begin{gathered} L(u_m(\mathbf{x})-v_m(\mathbf{x}))-\sum _{n=1}^N{s}_{mn}(u_n(\mathbf{x})-v_n(\mathbf{x}))-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}(u_m(\mathbf{x})-v_m(\mathbf{x})) \\ +\sum _{n=1}^N\frac{{\Psi }_n(0)}{f(\mathbf{x})}[u_n(\mathbf{x})(u_m(\mathbf{x})-v_m(\mathbf{x}))+v_m(\mathbf{x})(u_n(\mathbf{x})-v_n(\mathbf{x}))]=0. \end{gathered}$$

$$\left\{\begin{array}[]{rcll}\displaystyle\sum_{m=1}^{N}|{\rm div}(A\nabla h_{m})|^{2}&\leq&C_{1}\Big{(}\displaystyle\sum_{m=1}^{N}|h_{m}|^{2}+\sum_{m=1}^{N}|\nabla h_{m}|^{2}\Big{)}&\mbox{in }\Omega,\\ h_{m}&=&0&\mbox{on }\partial\Omega,\\ \partial_{\nu}h_{m}&=&0&\mbox{on }\partial\Omega\end{array}\right.$$

$$\sum _{m=1}^N{\int }_{\Omega }|\mathrm{div}(A\nabla h_m){|}^2\exp [2\lambda (x_d+b{)}^p]d\mathbf{x}\ge C\lambda \sum _{m=1}^N{\int }_{\Omega }\left[|\nabla h_m{|}^2+{\lambda }^2|h_m{|}^2\right]\exp \left[2\lambda {\left(x_d+b\right)}^p\right]d\mathbf{x}.$$

$$C_{1}\int_{\Omega}\Big{(}\displaystyle\sum_{m=1}^{N}|h_{m}|^{2}+\sum_{m=1}^{N}|\nabla h_{m}|^{2}\Big{)}\exp\left[2\lambda\left(x_{d}+b\right)^{p}\right]d{\bf x}\\ \geq C\lambda\sum_{m=1}^{N}\int_{\Omega}\left[|\nabla h_{m}|^{2}+\lambda^{2}|h_{m}|^{2}\right]\exp\left[2\lambda\left(x_{d}+b\right)^{p}\right]d\mathbf{x}.$$

$$L{v}_m^{(0)}(\mathbf{x})-\sum _{n=1}^N{s}_{mn}{v}_n^{(0)}(\mathbf{x})-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}{v}_m^{(0)}(\mathbf{x})=0\text{mbox}forall\mathbf{x}\in \Omega$$

$$v_m^{(0)}(\mathbf{x})=F_m(\mathbf{x}),{\partial }_{\nu }{v}_m^{(0)}(\mathbf{x})=G_m(\mathbf{x})\text{mbox}forall\mathbf{x}\in \partial \Omega$$

$$L{v}_m^{(p+1)}(\mathbf{x})-\sum _{n=1}^N{s}_{mn}{v}_n^{(p+1)}(\mathbf{x})-\frac{Lf(\mathbf{x})}{f(\mathbf{x})}{v}_m^{(p+1)}(\mathbf{x})+\sum _{n=1}^N\frac{{\Psi }_n(0)}{f(\mathbf{x})}{v}_n^{(p+1)}(\mathbf{x})v_m^{(p)}(\mathbf{x})=0\text{mbox}forall\mathbf{x}\in \Omega$$

Full text

Keywords: coefficient inverse problem, parabolic equations, approximation, Fourier coefficients

AMS Classification 35R30, 35K20

1 Introduction

Let $\Omega$ be a cube $(-R,R)^{d}\subset\mathbb{R}^{d}$, $d\geq 2$, where $R$ is a positive number. Introduce a $d\times d$ matrix valued function $A$ with entries in the class $C^{\infty}(\mathbb{R}^{d},\mathbb{R}^{d\times d})$. Assume throughout the paper that

$A$ is symmetric; i.e, $A^{T}=A$, 2. 2.

the matrix $A$ is uniformly elliptic; i.e., there exists a positive number $\mu$ such that

[TABLE] 3. 3.

for all ${\bf x}\in\mathbb{R}^{d}\setminus\Omega$, $A({\bf x})=\textrm{Id}$ where Id is the identity matrix.

Let ${\bf b}$ be a $d$-dimensional vector valued function in $C^{\infty}(\mathbb{R}^{d},\mathbb{R}^{d})$. Define the operator

[TABLE]

Consider the solution $u$ to the following initial value problem for the following parabolic equation

[TABLE]

Here, $f$ is a smooth function defined on $\mathbb{R}^{d}$. We refer the reader to [1, Chapter 7] and [2] for the well-posedness of (1.2). The regularity of the function $u({\bf x},t)$ can be found in those books. The second order term ${\rm div}(A({\bf x})\nabla u({\bf x},t))$ describes the diffusion, the first order term ${\bf b}({\bf x})\cdot\nabla u({\bf x},t)$ describes the transport and the zeroth order term $c({\bf x})u({\bf x},t)$ describes creation or depletion. In this paper, we numerically solve the problem of reconstructing the coefficient $c({\bf x})$ of the creation or depletion term.

Roughly speaking, the creation or depletion term refers to the ability of produce or destroy, respectively, photons. For e.g., in a chemical reaction, in which $u({\bf x},t)$ is the concentration of unstable gas, the rate of decomposition is proportional to the concentration. This leads to the presence of $c({\bf x})u({\bf x},t)$ in (1.2). We would like to refer the reader to [3, Chapter 6] for details about this term. In the current paper, we propose a method to solve the following highly nonlinear and severely ill-posed coefficient inverse problem.

Problem 1.1** (Coefficient Inverse Problem (CIP)).**

Assume that $f({\bf x})\not=0$ for all ${\bf x}\in\Omega$. Given a time $T>0$ and the lateral Neumann data

[TABLE]

for all ${\bf x}\in\partial\Omega$ and $t\in[0,T],$ determine the coefficient $c({\bf x})$, ${\bf x}\in\Omega$. Here $\nu({\bf x})$ is the outward normal vector of $\partial\Omega$ at a point ${\bf x}$.

Problem 1.1 has uncountable practical applications. In fact, suppose that interior points of a medium are not accessible. In this case, by measuring some boundary information of the function $u$, which are the heat and the heat flux in this paper, for a certain period of time and by solving Problem 1.1, one can determine the coefficient $c({\bf x})$ of the governing equation in (1.2), which enables us to inspect that medium without destructing it. We recall here a specific example in bioheat transfer. In this field, the coefficient $c({\bf x})$ represents the blood perfusion. The knowledge of this coefficient plays a crucial role in calculating the temperature of the blood flowing through the tissue, see [4]. The uniqueness of Problem 1.1 is still open and is studied in an approximation context of this paper. One can find the uniqueness of other versions of Problem 1.1 in [5, 6, 7] when some internal data are assumed to be known. When the Dirichlet to Neuman map is given, the reader can find the uniqueness in [8]. The uniqueness of Problem 1.1 is an assumption in this paper. Another related problem is the inverse problem of recovery other coefficients, for e.g., the diffusion, or the initial condition from the measurement of the final time for parabolic equations. This problem is very important and interesting, see [9, 10, 11, 12, 13] for theoretical results and numerical methods. In this paper, we model the function $u$ defined on $\mathbb{R}^{d}\times[0,\infty)$. However, the reader will see in Sections 2 and 4, our analysis and algorithm are made inside $\Omega$. In other words, our method works when domain of the function $u$ and model (1.2) is restricted to $\Omega$.

Coefficient inverse problems for parabolic equations were studied intensively. Up to the knowledge of the author, the widely used method to solve this problem is the optimal control approach, see e.g., [14, 15, 4, 16, 17] and references therein. The authors of [14] applied the optimal control method involving a preconditioner to numerically compute the thermal conductivity with high quality. The main drawback of this method is the need of a good initial guess for the true solution which is not always available. On the other hand, we specially draw the reader’s attention to the convexification method, see [18, 19], which can overcome the difficulty about the availability of the initial guess. In those papers [18, 19], the authors introduce a convex functional whose minimizer yields the solution of the problem under consideration, by combining the quasi-reversibility method and the Carleman weight functions. One dimensional numerical examples, illustrating the role of Carleman weight functions in convexifying the cost functionals, are presented in [18]. It is valuable to numerically test this convexification method in higher dimensions. We also cite to [20] for another method to solve Problem 1.1 by repeatedly solving its linearization. In the current paper, we propose a novel method in which no advanced knowledge about the true coefficient is required.

Our method to solve Problem 1.1 consists of two stages. In the first stage, we eliminate the function $c({\bf x})$ from (1.2). The resulting equation obtained in this stage is not a standard equation. A numerical method to solve it is not available yet. In the second stage, we approximate that nonstandard equation as a coupled system of elliptic partial differential equations. This system is derived based on a truncation of the Fourier series, with respect to a special basis originally introduced in [21]. We apply a predictor-corrector procedure, in which the first approximation of the true solution is computed without any of its advance knowledge. The solution of Problem 1.1 follows.

Two important steps in our method require us to find vector valued functions satisfying a system of elliptic partial differential equations and both Dirichlet and Neumann boundary conditions. We employ the quasi-reversibility method, so-called a global minimization, for this purpose and we also prove the convergence of the quasi-reversibility method in our context, using a new Carleman estimate in [22]. The quasi-reversibility method was first introduced by Lattès and Lions in [23] for numerical solutions of ill-posed problems for partial differential equations. It has been studied intensively since then, see e.g., [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 22]. A survey on this method can be found in [34].

In Section 2, we derive the system mentioned above. In Section 3, we study the uniqueness of this system. In Section 4, we propose a numerical method to solve that system. Also in Section 4, we study the quasi-reversibility method that can be applied in our context. In Section 5, we describe the implementation using the finite difference method. In Section 6, we present some numerical results. Section 7 is for the concluding remarks.

2 A nonlinear coupled system of elliptic equations

In this paper, the initial condition of $f$ is assumed to be strictly positive in $\overline{\Omega}$ and in $H^{2+\beta}(\mathbb{R}^{d})$ for some $\beta>0$ and has compact support in $\mathbb{R}^{d}$. Then, since all coefficients of the operator $L$ are in class $C^{\infty}$, (1.2) has a unique solution with $|u({\bf x},t)|\leq M$ and $u\in H^{2+\beta,1+\beta/2}(\mathbb{R}^{d}\times[0,T])$ for some constant $M>0$. These unique solvability and regularity properties can be obtained by applying Theorem 6.1 in [2, Chapter 5, §6] and Theorem 2.1 in [2, Chapter 5, §2]. We impose the condition that $\beta$ is large such that the function $u$ has second derivative with respect to time. Further more, if we impose a stronger condition that $f$ belongs to the class $C_{0}^{\infty}(\mathbb{R}^{d})$, then applying the induction arguments in the third paragraph of [2, page 456], we see that the function is $k$ times differentiable for any $k>0$.

From now on, we denote by $\Omega_{T}$ the set $\Omega\times[0,T].$ Define the function

[TABLE]

It follows from (1.2) that

[TABLE]

On the other hand, for all ${\bf x}\in\Omega,$

[TABLE]

Therefore,

[TABLE]

Plugging (2.3) into (2.2), we obtain the following equation

[TABLE]

Remark 2.1**.**

Solving the nonlinear equation (2.4) is challenging due to the presence of the initial condition $v({\bf x},0)$. A theoretical result to solve it is not yet available. We employ the technique of truncating the Fourier series, see [21], to solve (2.4).

Recall a special orthonormal basis of $L^{2}(0,T)$ originally introduced by Klibanov [21] in 2017. This basis plays a crucial role in deriving an approximate model whose solution will be used to directly compute the solution of Problem 1.1. For each $n\geq 1$, define the function $\phi_{n}(t)=(t-T/2)^{n-1}\exp(t-T/2)$. It is well-known that the set $\{\phi_{n}\}_{n=1}^{\infty}$ is complete in $L^{2}(0,T)$. Employing the Gram-Schmidt orthonormalization process on this set, we obtain an orthonormal basis of $L^{2}(0,T)$. We denote this basis by $\{\Psi_{n}\}_{n=1}^{\infty}.$ We have the following proposition.

Proposition 2.1** (see [21]).**

The basis $\{\Psi_{n}\}_{n=1}^{\infty}$ satisfies the following properties:

$\Psi_{n}$* is not identically zero for all $n\geq 1$,* 2. 2.

For all $m,n\geq 1$

[TABLE]

As a result, for all integer $N>0$, the matrix $S=(s_{mn})_{m,n=1}^{N}$, is invertible.

Recall the Fourier coefficients of the function $v({\bf x},t)$

[TABLE]

We have

[TABLE]

Fix a number $N>0$. We approximate the function $v({\bf x},t)$ by the partial sum

[TABLE]

In this approximation context,

[TABLE]

and

[TABLE]

Plugging (2.6), (2.7) and (2.8) into (2.4), we have

[TABLE]

for all $({\bf x},t)\in\Omega_{T}.$ For each $m$ in $\{1,2,\dots,N\}$, multiply both sides of the equation above by $\Psi_{m}(t)$ and then integrate the resulting equation with respect to $t$. Noting that

[TABLE]

we have for each $m\in\{1,\dots,N\}$,

[TABLE]

for all ${\bf x}\in\Omega.$ On the other hand, for each $m\in\{1,\dots,N\}$, the function $v_{m}({\bf x})$ satisfies the following constraints for all ${\bf x}\in\partial\Omega$

[TABLE]

Remark 2.2** (Data).**

From now on, we consider $F_{m}$ and $G_{m}$ as the Cauchy indirect data. These data are computed directly from the time dependent data via (2.10).

Remark 2.3** (Noise).**

Let $F^{*}$ and $G^{*}$ be the data without noise. Let ${\rm rand}$ be the function of uniformly distributed random numbers in the range $[-1,1]$. For $\delta>0$, define

[TABLE]

The corresponding noisy indirect data are denoted by $F^{\delta}_{m}$ and $G^{\delta}_{m}$ respectively. There is a difficulty in computing $F^{\delta}_{m}$ and $G^{\delta}_{m}$ via (2.10). Computing the derivatives $F_{t}$ and $G_{t}$ of the data for (2.10) is unstable. We use the well-known Tikhonov regularization technique to compute them. In this paper, we test our numerical method for simulated data with noise level $\delta=5\%$.

In summary, we have proved the following proposition.

Proposition 2.2**.**

Fix $N>0$. Assume that the function $v({\bf x},t)$, for $({\bf x},t)\in\Omega_{T}$, can be well-approximated by the expression in (2.6) with the function $v_{n}({\bf x})$, $n\in\{1,2,\dots,N\}$, given in (2.5). Then the Fourier coefficients $v_{n}$ satisfies the over-determined system of partial differential equations (2.9)–(2.10).

3 The uniqueness of Problem 1.1 in the approximation context

Proposition 2.2 suggests a method to numerically solve Problem 1.1. We solve the nonlinear system (2.9)–(2.10) for a vector $(v_{1},\dots,v_{N})$ and compute the coefficient $c$ via (2.6) and then (2.3). Hence, the uniqueness of the solution to the nonlinear system (2.9)–(2.10) implies the unique reconstruction of Problem 1.1 assuming the approximation (2.6). In this section, we establish the uniqueness for (2.9)–(2.10). It is sufficient to study this uniqueness for the bounded solution. In fact, due to the conditions imposed on the source function $f$ in the first paragraph of Section 2, namely $f\in C_{0}^{\infty}(\mathbb{R}^{d})$, the true function $u({\bf x},t)$ is $k$ times differentiable for any $k>0$. Thus, the function $v({\bf x},t)=u_{t}({\bf x},t)$ is bounded. It follows from (2.5) that the true solution to (2.9)–(2.10) is in $L^{\infty}(\Omega)$.

In order to establish the uniqueness of the bounded solution to (2.9)–(2.10), we need the following Carleman estimate.

Lemma 3.1** (Carleman estimate).**

Let the number $b>R$. Then there exist numbers $p_{0}\geq 1$ and $\lambda_{0}\geq 1$ depending only on $\mu$, $b$, $d$, $R$, $\|A\|_{L^{\infty}(\Omega)^{d\times d}}$ such that the following Carleman estimate holds:

[TABLE]

for all $\lambda\geq\lambda_{0},$ $p\geq p_{0}$ and $u\in H^{2}(\Omega)$ with $u=\partial_{\nu}0=0$ on $\partial\Omega$. Here, the constant $C$ depends only on $\mu$, $b$, $d$, $R$ and $\|A\|_{L^{\infty}(\Omega)^{d\times d}}$.

Lemma 3.1 is a direct consequence of [22, Theorem 4.1]. We do not repeat the proof in this paper. This Lemma plays a crucial role in the proof of Theorem 3.1 below. Moreover, it will be applied to prove the convergence of the quasi-reversibility method, see Theorem 4.9.

We have the theorem.

Theorem 3.1**.**

The nonlinear system (2.9)–(2.10) has at most one solution in $H^{2}(\Omega)\cap L^{\infty}(\Omega)$.

Proof.

Let $U=(u_{1},\dots,u_{N})$ and $V=(v_{1},\dots,v_{N})$ be two solutions in $H^{2}(\Omega)\cap L^{\infty}(\Omega)$ to (2.9)–(2.10). It follows from (2.9) that for all ${\bf x}\in\Omega$ and $m\in\{1,\dots,N\}$,

[TABLE]

and

[TABLE]

It follows from (3.2) and (3.3) that for all ${\bf x}\in\Omega$ and $m\in\{1,\dots,N\}$,

[TABLE]

Let $h=(h_{1},\dots,h_{N})=U-V$. Since $U$ and $V$ are in $L^{\infty}(\Omega),$ we can find a number $M$ such that $|u_{m}|\leq M$ and $|v_{m}|\leq M$ for all $m\in\{1,\dots,N\}.$ It follows from (1.1), (2.10) and (3.4) that

[TABLE]

for some constant $C_{1}$ depending only on $N$, $M$, $(s_{mn})_{m,n=1}^{N}$, $(\Psi_{m})_{m=1}^{N}$, $f$, and ${\bf b}$. Let $b>0$, $\lambda>\lambda_{0}$ and $p>p_{0}$ where $\lambda_{0}$ and $p_{0}$ are in Lemma 3.1. Applying the Carleman estimate in Lemma 3.1 for each function $h_{m}$, $m\in{1,\dots,N}$, we have

[TABLE]

It follows from (3.5) and (3.6) that

[TABLE]

Choosing $\lambda$ sufficiently large, we obtain $h_{m}=0$ for all $m\in\{1,\dots,N\}.$ ∎

Remark 3.1**.**

The uniqueness of the solution $(v_{1},\dots,v_{N})$ to the nonlinear system (2.9)–(2.10) implies the unique reconstruction of the solution to Problem 1.1. This can be seen via the reconstruction formulas (2.6) and (2.3). This uniqueness holds true only in the approximation context (2.6) while the uniqueness for the true model in the time domain is extremely challenging, which is out of the scope of this paper.

4 The method to solve Problem 1.1

Solving Problem 1.1 is reduced to solving the system of nonlinear partial differential equations (2.9)–(2.10).

4.1 An iterative process

We propose the following iterative method, in which a predictor-corrector procedure is applied. The first predictor, named as $V^{(0)}$, is set to be the solution of the linear system obtained by removing from (2.9) the nonlinear term. More precisely, we set $V^{(0)}=(v_{1}^{(0)},\dots,v_{N}^{0})^{T}$ as the solution of

[TABLE]

and

[TABLE]

for $m\in\{1,\dots,N\}.$ Next, by induction, assume that $V^{(p)}$ is known for some positive integer $p$, we find $V^{(p+1)}$ by solving the equation obtained from (2.9) by replacing $v_{m}$ in the nonlinear term by its approximation $v_{m}^{(p)}.$ That means, $V^{(p+1)}=(v_{1}^{(p+1)},\dots,v_{N}^{(p+1)})^{T}$ is set to be the solution of

[TABLE]

and

[TABLE]

for each $m\in\{1,\dots,N\}.$

Due to the presence of the latteral Cauchy data, both problems (4.1)–(4.2) and (4.3)–(4.4) are over-determined. In the case when the data is noisy, they might not have a solution. We, therefore, employ the quasi-reversibility method to solve them. In Proposition 4.1, we show the existence of an “approximation” of the true solution in the case when the data has no noise. In Theorem 4.9, we prove the convergence of this approximation as the noise tends to [math].

4.2 The quasi-reversibility method

We next recall the quasi-reversibility method to solve systems of partial differential equations with Cauchy boundary data. The two systems of elliptic partial differential equations (4.1)–(4.2) and (4.3)–(4.4) are over-determined due to both Dirichlet and Neumann boundary conditions imposed. We use the quasi-reversibility method to solve them. A general form of the problem (4.1)–(4.2) and the problem (4.3)–(4.4) is given by

[TABLE]

where $A$ is introduced in Section 1 and $B$ is a $N\times N$ matrix valued function in $L^{\infty}(\Omega).$ We have the following proposition whose proof closely follows that of Theorem 3.1 in [22].

Proposition 4.1**.**

Fix $\epsilon>0$. Then, the functional

[TABLE]

has a unique minimizer on $H^{2}(\Omega).$ This minimizer $V_{\epsilon}$ is called the regularized solution of (4.5).

In summary, we propose Algorithm 1 to solve Problem 1.1 via solving (2.9)–(2.10) by the quasi-reversibility method.

Proposition 4.2** (The uniqueness of (4.5)).**

The system (4.5) has at most one solution in $H^{2}(\Omega)$.

The proof of this proposition is very similar to that of Theorem 3.1 by using the Carleman estimate. We do not repeat the proof here.

Theorem 4.1** (The convergence of the quasi-reversibility method for (4.5)).**

Assume that there uniquely exists a true solution $V^{*}$ to (4.5) with the boundary data $\mathcal{F}$ and $\mathcal{G}$ replaced by the corresponding noiseless ones, denoted by $\mathcal{F}^{*}$ and $\mathcal{G}^{*}$ respectively. Let $\mathcal{F}^{\delta}$ and $\mathcal{G}^{\delta}$ be the corresponding noisy data for some $\delta>0$. Assume that there exists an “error” vector valued function $\mathcal{E}$ such that

[TABLE]

and assume that

[TABLE]

Then, $V_{\epsilon}^{\delta}$, the minimizer of $J_{\epsilon}$ with $\mathcal{F}$ and $\mathcal{G}$ replaced by $\mathcal{F}^{\delta}$ and $\mathcal{G}^{\delta}$ respectively, satisfies the estimate

[TABLE]

Remark 4.1**.**

The convergence for the quasi-reversibility method guaranteed by Theorem 4.9 is similar to that in [22, Theorem 5.1]. The main difference of two results is in the objective functional is minimized subject to some boundary constraints while in the current paper, such constraints are relaxed by adding the two boundary integrals in (4.6).

Remark 4.2** (The existence of the error function $\mathcal{E}$ in the statement of Theorem 4.9).**

By employing the main result in [26], we can prove the existence of the vector valued function $\mathcal{E}$ that satisfies (4.7) if we impose the reasonable conditions that $\Omega$ is in the class $C^{1,1}$ and that $(\mathcal{F}^{\delta}-\mathcal{F}^{*},\mathcal{G}^{\delta}-\mathcal{G}^{*})$ belongs to $H^{3/2}(\partial\Omega)$.

Proof of Theorem 4.9.

Since $V_{\epsilon}^{\delta}$ is the regularized solution to (4.5), it is the minimizer of $J_{\epsilon}$, defined in (4.6). Hence, for all $\phi\in H^{2}(\Omega)^{N}$, we have

[TABLE]

for all $\phi\in H^{2}(\Omega)^{N}.$ On the other hand, since $V^{*}$ is the true solution to (4.5)

[TABLE]

for all $\phi\in H^{2}(\Omega)^{N}.$ Taking the difference of (4.10) and (4.11), we have

[TABLE]

where $W=V^{\delta}_{\epsilon}-V^{*}$ for all $\phi\in H^{2}(\Omega)^{N}.$ Using

[TABLE]

as a test function in (4.12) and using (4.7), we have

[TABLE]

Applying the inequality $|\langle u,v\rangle|\leq 1/2(\|u\|^{2}+\|v\|^{2})$, (4.8) and the trace theory, we have

[TABLE]

Here, $C$ is a generic constant that might change from estimate to estimate. Choose $b>R$, $\lambda>\lambda_{0}$, $p\geq p_{0}$ where $b,$ $\lambda_{0}$ and $p_{0}$ are as in Lemma 3.1. It is not hard to verify that the function $\phi$ satisfies the homogenous boundary conditions $\phi=\partial_{\nu}\phi=0$ on $\partial_{\Omega}\times[0,T]$. Using the Carleman estimate in Lemma 3.1, we can bound the left hand side of (4.14) as follows

[TABLE]

Choosing $\lambda$ sufficiently large, since $B\in L^{\infty}(\Omega),$ we have

[TABLE]

This, together with (4.13) and (4.14), implies (4.9). The theorem is proved. ∎

Corollary 4.1**.**

Theorem 4.9 implies that the regularized solution to (4.5) obtained by quasi-reversibility method is close to the true solution. In fact, if $\epsilon<O(\delta^{2})$, it follows from (4.9) that

[TABLE]

The rate of convergence is Lipschitz.

Theorem 4.9 guarantees that Steps 2 and 4 in Algorithm 1 below provide good approximations of the sequence $\{c_{p}^{\delta}\}_{p=1}^{\infty}$ in comparison to the sequence $\{c_{p}^{*}\}_{p=1}^{\infty}$ with the Lipschitz rate provided that $\epsilon=O(\delta^{2})$ as $\delta\to 0^{+}$. If the sequence $\{c_{p}^{*}\}_{p=1}^{\infty}$ converges to the solution of Problem 1.1, Algorithm 1 yields a numerical procedure to solve it. This convergence is verified numerically in Section 6.

4.3 The procedure to solve the coefficient inverse problem for parabolic equations

By Proposition 2.2, the strategy to solve (2.9)–(2.10) described in Section 4.1 and the convergence of the quasi-reversibility method, see Theorem 4.9, we propose Algorithm 1 to reconstruct the coefficient $c({\bf x})$ for ${\bf x}\in\Omega.$

Remark 4.3**.**

Unlike the widely used least squares method to solve ill-posed inverse problems, we do not require a good initial guess for the true coefficient $c({\bf x})$. Our first approximation is computed in Steps 2 and 3 of Algorithm 1. It is shown in Section 6 that the functions $c^{(0)}$ are acceptable in Tests 1, 3 and 4. In contrast, our computed $c^{(0)}$ is poor in Test 2. However, the error is automatically corrected when we find $c^{(1)}$ in Steps 4 and 5.

5 The implementation using the finite difference method

We test our method in the simple case when $d=2$, $\Omega=(-1,1)^{2}$, $L=\Delta$ is the Laplacian and $f({\bf x})=100$.

5.1 The forward problem

To generate the simulated data, we solve the forward problem of Problem 1.1. That means, given $c({\bf x})$ (see each test below for the definition of $c({\bf x})$), we need to compute the solution $u({\bf x},t)$ to (1.2) on the whole plane $\mathbb{R}^{2}$. Instead of doing so, we solve an analog of (1.2) on a domain $\Omega_{1}=({-}R_{1},R_{1})^{2}$ where $R_{1}=3>R=1.$

This domain approximation does not effect our analysis because the formulation of the parabolic equation in (1.2) is used only inside the domain $\Omega$ when we study the inverse problem.

In other words, we solve the equation

[TABLE]

Here, we choose the time-independent Dirichlet boundary data for the simplicity. In this paper, we solve problem (5.1) by implicit method using finite differences by the backward Euler scheme. In the finite difference scheme, we find the function $u({\bf x},t)$ on the grid of points

[TABLE]

where $N_{1}$ and $N_{t}$ are two large integers, $d_{{\bf x}_{1}}=2R_{1}/(N_{1}-1)$ and $d_{t}=T/(N_{t}-1)$. In our computational program, $N_{1}$ is set to be 240 and $N_{t}=100.$ Having the function $u({\bf x},t)$ for all ${\bf x}\in\Omega_{1}$ and $t\in[0,T]$ in hand, we can directly extract the data $F({\bf x},t)=u({\bf x},t)$ and $G({\bf x},t)=\partial_{\nu}u({\bf x},t)$ on $\partial\Omega\times[0,T]$.

5.2 The inverse problem

In this section, we present how to implement Algorithm 1 in the finite difference scheme. Similarly to the previous section, we define and then compute the function $c$ on a uniform grid of points

[TABLE]

where $N_{\bf x}$ and $N_{t}$ are two large integers, $d_{\bf x}=2R/(N_{\bf x}-1)$ and $d_{t}=T/(N_{t}-1)$. In all numerical tests in Section 6, we take $N_{\bf x}=80$ and $N_{t}=100.$

We next present each step of Algorithm 1.

Step 1. In this step, to choose “truncation” number $N$. To do so, we take a “reference” function $v$ in one of the examples in Section 6 and then compute the absolute difference

[TABLE]

We observe that the larger $N$, the smaller $\|e_{N}\|_{L^{\infty}(\Omega\times[0,T])}$. We examine the function $e_{N}$ when $N=5,$ $N=10$ and $N=25$, see Figure 1. It is evident from Figure 1c that when $N=25$, $\|e_{N}\|_{L^{\infty}}$ is sufficiently small, about $4{\times}10^{-3}$.

As a result, we choose $N=25$. We use this choice of $N$ for all numerical tests. We observe that using higher $N$ does not improve the quality of the reconstructed coefficient $c({\bf x})$. Also in Step 1 of Algorithm 1, we choose the regularized parameter $\epsilon=10^{-9}.$ This choice of $\epsilon$ is based on a trial-error process in the case when the given data is noisy.

Step 2. Compute the vector valued function $V^{(0)}$, which is set to be the minimizer of the functional, due to the quasi-reversibility method,

[TABLE]

Here for each $m\in\{1,\dots,N\},$ $F_{m}=F_{m}^{\delta}$ and $G_{m}=G_{m}^{\delta}$ are computed in (2.10) with $F^{\delta}$ and $G^{\delta}$ replacing $F$ and $G$ respectively with $\delta=5\%$.

Here, we replace the term $\|V\|_{H^{2}(\Omega)}^{2}$ in (4.6) by the term $\|V\|_{H^{1}(\Omega)}^{2}.$ This is because the $H^{1}(\Omega)-$norm is easier to work with computationally than the $H^{2}(\Omega)-$norm. On the other hand, we have not observed any instabilities probably because the number $80\times 80$ of grid points we use is not too large and all norms in finite dimensional spaces are equivalent. We now identify $\{v_{m}(x_{i},y_{j}):1\leq i,j,N_{\bf x},1\leq m\leq N\}$ by the $N_{\bf x}^{2}N$ dimensional vector $\mathfrak{v}$ whose $\mathfrak{i}^{\rm th}$ entry is given by

[TABLE]

Here, $(i,j,m)$ is such that

[TABLE]

Then, by approximating all differential operators in the right hand side of (5.3) with their corresponding finite difference versions, we have

[TABLE]

where the matrices $\mathfrak{L}$, $\mathfrak{D}_{1}$, $\mathfrak{D}_{2}$, $D_{x}$ and $D_{y}$ and the vectors $\mathfrak{F}$ and $\mathfrak{G}$ are described below. The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $\mathfrak{L}$ is given by

$\mathfrak{L}_{\mathfrak{i}\mathfrak{j}}=-4/d_{\bf x}^{2}-s_{mm}-\Delta f(x_{i},y_{j})/f(x_{i},y_{j})$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$; 2. 2.

$\mathfrak{L}_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}^{2}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i\pm 1-1)N_{\bf x}N+t(j\pm 1-1)N+m$; 3. 3.

$\mathfrak{L}_{\mathfrak{i}\mathfrak{j}}=-s_{mn}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+n$, $n\not=m;$ 4. 4.

all other entries of $\mathfrak{L}$ are [math];

for all $2\leq i,j\leq N_{\bf x}-1,$ $1\leq m\leq N.$ The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $\mathfrak{D}_{1}$ is given by

$(\mathfrak{D}_{1})_{\mathfrak{i}\mathfrak{j}}=1$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $i\in\{1,N_{\bf x}\}$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 2. 2.

$(\mathfrak{D}_{1})_{\mathfrak{i}\mathfrak{j}}=1$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $1\leq i\leq N_{\bf x},$ $j\in\{1,N_{\bf x}\}$, $1\leq m\leq N;$ 3. 3.

all other entries of $\mathfrak{D}_{1}$ are [math].

The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $\mathfrak{D}_{2}$ is given by

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $i\in\{1,N_{\bf x}\}$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 2. 2.

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i+1-1)N_{\bf x}N+(j-1)N+m$ for $i=1$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 3. 3.

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i-1-1)N_{\bf x}N+(j-1)N+m$ for $i=N_{\bf x}$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 4. 4.

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $i1\leq i\leq N_{\bf x}$, $j\in\{1,N_{\bf x}\}$, $1\leq m\leq N;$ 5. 5.

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i-1)N_{\bf x}N+(j+1-1)N+m$ for $1\leq i\leq N_{\bf x}$, $j=1$, $1\leq m\leq N;$ 6. 6.

$(\mathfrak{D}_{2})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i-1)N_{\bf x}N+(j-1-1)N+m$ for $1\leq i\leq N_{\bf x}$, $j=N_{\bf x}$, $1\leq m\leq N;$ 7. 7.

all other entries of $\mathfrak{D}_{2}$ are [math].

The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $D_{x}$ is given by

$(D_{x})_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $1\leq j\leq N_{\bf x}-1$, $1\leq j\leq N_{\bf x}-1$, $1\leq m\leq N;$ 2. 2.

$(D_{x})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i+1-1)N_{\bf x}N+(j-1)N+m$ for $i=1$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 3. 3.

all other entries of $D_{x}$ are [math].

The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $D_{y}$ is given by

$(D_{y})_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$ for $1\leq j\leq N_{\bf x}-1$, $1\leq j\leq N_{\bf x}-1$, $1\leq m\leq N;$ 2. 2.

$(D_{y})_{\mathfrak{i}\mathfrak{j}}=-1/d_{\bf x}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j+1-1)N+m$ and $\mathfrak{j}=(i+1-1)N_{\bf x}N+(j-1)N+m$ for $i=1$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 3. 3.

all other entries of $D_{y}$ are [math].

The vector $\mathfrak{F}$ is defined as

$\mathfrak{F}_{\mathfrak{i}}=F_{m}(x_{i},y_{j})$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ for $i\in\{1,N_{\bf x}\}$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 2. 2.

all other entries of $\mathfrak{F}$ are [math].

The vector $\mathfrak{G}$ is defined as

$\mathfrak{G}_{\mathfrak{i}}=G_{m}(x_{i},y_{j})$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ for $i\in\{1,N_{\bf x}\}$, $1\leq j\leq N_{\bf x}$, $1\leq m\leq N;$ 2. 2.

all other entries of $\mathfrak{G}$ are [math].

Since $\epsilon$ is small (in our computational program $\epsilon=10^{-9}$), to find the minimizer $V^{(0)}$ of the finite difference version of $J_{\epsilon}$ in (5.6), we solve the linear system

[TABLE]

Having $\mathfrak{v}$ in hand, we can compute $V^{(0)}=(v_{1}^{0},\dots,v_{N}^{(0)})^{\rm T}$ using (5.4) and (5.5).

Step 4. The implementation for this step is similar to that for Step 2, namely, we minimize

[TABLE]

To this end, we identify the vector valued function $V$ by the vector $\mathfrak{v}$ as in (5.4) and (5.5) and then solve the linear system

[TABLE]

Here, the matrices $\mathfrak{D}_{1}$, $\mathfrak{D}_{2}$, $D_{x}$ and $D_{y}$ and the vectors $\mathfrak{F}$ and $\mathfrak{G}$ are defined in the implementation section for Step 2. The $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$ matrix $\mathcal{L}$ is given by

$\mathcal{L}_{\mathfrak{i}\mathfrak{j}}=-4/d_{\bf x}^{2}-s_{mm}-\Delta f(x_{i},y_{j})/f(x_{i},y_{j})+v_{m}^{(p)}(x_{i},y_{j})\Psi_{m}(0)/f(x_{i},y_{j})$ if $\mathfrak{i}=\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+m$; 2. 2.

$\mathcal{L}_{\mathfrak{i}\mathfrak{j}}=1/d_{\bf x}^{2}$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i\pm 1-1)N_{\bf x}N+(j\pm 1-1)N+m$; 3. 3.

$\mathcal{L}_{\mathfrak{i}\mathfrak{j}}=-s_{mn}+v_{n}^{(p)}(x_{i},y_{j})\Psi_{n}(0)/f(x_{i},y_{j})$ if $\mathfrak{i}=(i-1)N_{\bf x}N+(j-1)N+m$ and $\mathfrak{j}=(i-1)N_{\bf x}N+(j-1)N+n$, $n\not=m;$ 4. 4.

all other entries of $\mathcal{L}$ are [math];

for all $2\leq i,j\leq N_{\bf x}-1,$ $1\leq m\leq N.$ Having $\mathfrak{v}$ in hand, we can compute $V^{(p+1)}=(v_{1}^{(p+1)},\dots,v_{N}^{(p+1)})^{T}$ using (5.4) and (5.5).

Steps 3, 5 and 6. The implementation of these steps is straight forward.

Remark 5.1**.**

All matrices above are of the large size $N_{\bf x}^{2}N\times N_{\bf x}^{2}N$, which might cause some inefficiency in computations. However, since most of their entries are zeros, we can treat those matrices as sparse ones to overcome this difficulty. The linear algebra package for sparse matrices are already built in Matlab.

In the next section, we show some numerical results.

6 Numerical examples

The numerical results presented below are computed from the knowledge of $F_{m}({\bf x})$ and $G_{m}({\bf x})$, $m\in\{1,\dots,N\}$, on $\partial\Omega\times[0,0.3]$ including 10% of noise, where $F$ and $G$ are the boundary data in Remark 2.3. The number of truncation $N$ is $25$. The regularization parameter is $\epsilon=10^{-9}.$ The computational program is implemented by the finite difference method.

Test 1. The true function $c_{\rm true}$ has a smooth inclusion

[TABLE]

The numerical results for this case are displayed in Figure 2. One can observe in Figures 2b–2e that the circular shape and location of the inclusion can be succesfully detected. The true maximal value of the function $c_{\rm true}$ is 20. The reconstructed maximal value of the function $c_{\rm comp}=c^{(10)}$ is $19.07$. 2. 2.

Test 2. We test the case when the function $c_{\rm true}$ is a step function with two rectangular inclusions. This example is interesting since $c_{\rm true}$ is not smooth and the gap at the boundaries of the inclusions is high. The function $c_{\rm true}$ is given by

[TABLE]

The numerical results for this case are displayed in Figure 3. One can observe in Figures 3c–3e that the reconstructed rectangular shape and location of the inclusion are satisfactory. The true maximal value of the function $c_{\rm true}$ is 10. The reconstructed maximal value of the function $c_{\rm comp}$ is 10.98. Similarly to the previous test, it is evident from Figure 3f that our method converges fast. 3. 3.

Test 3. We test the case of two circular inclusions. In this case, the function $c_{\rm comp}$ is a step function with high gap at the boundary of the inclusions. The function $c_{\rm true}$ is given by

[TABLE]

The numerical results for this case are displayed in Figure 4. One can observe in Figure 4b that the circular shape and can be successfully detected at the first step. The true maximal value of the function $c_{\rm true}$ at the lower inclusion is 8 and the reconstructed one is 8.90. The true maximal value of the function $c_{\rm true}$ at the upper inclusion is 5 and the reconstructed one is 5.24. Figure 4f shows the stability of our method. 4. 4.

Test 4. We test the case when the function $c_{\rm true}$ is allowed to be negative. In this case, the function $c_{\rm comp}$ is the letter $X$ with a half is positive and another half is negative. The function $c_{\rm true}$ is given by

[TABLE]

The numerical results for this case are displayed in Figure 5. The reconstructed image of the letter $X$ is acceptable. The true maximal positive value of the function $c_{\rm true}$ is 8 and the reconstructed one is 8.90. The true minimal negative value of the function $c_{\rm true}$ is -8 and the reconstructed one is $-7.93$. Again, Figure 5f shows the stability of our method.

Remark 6.1**.**

It is evident from Figures 2c–5c that Algorithm 1 is robust in the sense that it provides good reconstructed coefficient $c_{\rm comp}$ after a few iterations without any requirement of an initial guess. It is remarkable mentioning that, in all tests above, our method provides good numerical results without any advanced knowledge of the true coefficient $c_{\rm true}$. However, as seen in 2e–5e, there are some artifacts. These artifacts might be caused by cutting off the Fourier series of the function $v({\bf x},t)$ in (2.6).

7 Concluding remarks

In this paper, we introduced a new approach to numerically compute the creation or depletion coefficient of a general parabolic equation from lateral Cauchy data. Although this problem is highly nonlinear, we successfully compute this coefficient without requiring a good initial guess. In the first step, we derive an equation without the presence of the unknown coefficient. Then, we consider an approximation of the solution to this equation by truncating its Fourier series with respect to a special orthonormal basis of $L^{2}$. By this, we obtain a system of nonlinear elliptic equations. Numerically solving this system by an iterative procedure directly yields the desired coefficient. Two of the important strengths of this method are that, unlike the optimal control method, we do not require a good initial guess for the true solution and that our algorithm converges fast. However, the main drawback in this paper is the analysis of the convergence of the sequence obtained by the algorithm. A theorem that guarantees the efficiency of our method is not yet available. That means, our method is verified only in the numerical point of view. On the other hand, we strongly believe that our technique can be applied to compute the diffusion coefficients of parabolic equations. This serves as our near future work.

Acknowledgments

The author is grateful to Michael V. Klibanov for many fruitful discussions.