A numerical method for an inverse source problem for parabolic equations and its application to a coefficient inverse problem

TL;DR

This paper introduces a two-stage numerical method combining derivative equations and quasi-reversibility to solve inverse source and nonlinear coefficient inverse problems for parabolic equations, with demonstrated numerical results.

Contribution

It develops a novel two-stage numerical approach for inverse source problems and applies it iteratively to solve nonlinear coefficient inverse problems, advancing computational methods in this area.

Findings

Successful reconstruction of sources from observations.

Effective iterative solution for nonlinear coefficient inverse problems.

Numerical validation demonstrating method accuracy.

Abstract

Two main aims of this paper are to develop a numerical method to solve an inverse source problem for parabolic equations and apply it to solve a nonlinear coefficient inverse problem. The inverse source problem in this paper is the problem to reconstruct a source term from external observations. Our method to solve this inverse source problem consists of two stages. We first establish an equation of the derivative of the solution to the parabolic equation with respect to the time variable. Then, in the second stage, we solve this equation by the quasi-reversibility method. The inverse source problem considered in this paper is the linearization of a nonlinear coefficient inverse problem. Hence, iteratively solving the inverse source problem provides the numerical solution to that coefficient inverse problem. Numerical results for the inverse source problem under consideration and the corresponding nonlinear coefficient inverse problem are presented.

Figures

Tables

Table 1. Table 1: Test 2. The local extreme values of the functions p true subscript 𝑝 true p_{\rm true} and p comp subscript 𝑝 comp p_{\rm comp} at two inclusions. The relative error is denoted by err rel subscript err rel {\rm err}_{\rm rel} . Inclusion 1 is the one on the top right and inclusion 2 is the one on the bottom left.

Inclusion	noise level	extreme valuetrue	extreme valuecomp	$\left({\mathrm{err}}_{\mathrm{rel}}\right)$
1	0%	1	0.988	1.2%
2	0%	-2	-1.99	0.5%
1	5%	1	0.984	1.6%
2	5%	-2	-1.987	0.65%
1	10%	1	1.004	0.4%
2	10%	-2	-1.938	3.1%

Table 2. Table 2: Test 3. The true and reconstructed extreme values of inclusions. Inclusion 1 is the ellipse and inclusion 2 is the disk. The relative error is denoted by err rel subscript err rel {\rm err}_{\rm rel} .

Inclusion	Noise level	Extreme valuetrue	Extreme valuecomp	${\mathrm{err}}_{\mathrm{rel}}$
1	0%	2	2.151	7.5%
2	0%	-2	-2.133	6.7%
1	5%	2	2.211	10.6%
2	5%	-2	-2.117	5.9%
1	10%	2	2.312	15.6%
2	10%	-2	-2.29	14.5%%

Equations

$$\left\{\begin{array}[]{rcll}u_{t}({\bf x},t)&=&\mathcal{A}u({\bf x},t)+f({\bf x},t)p({\bf x})&{\bf x}\in\Omega,t>0,\\ u({\bf x},t)&=&0&{\bf x}\in\partial\Omega,t>0,\\ u({\bf x},0)&=&0&{\bf x}\in\Omega\end{array}\right.$$

$$G(\mathbf{x},t)={\partial }_{\nu }u(\mathbf{x},t)$$

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

$$v_t(\mathbf{x},t)=\mathcal{A}v(\mathbf{x},t)+f_t(\mathbf{x},t)p(\mathbf{x})$$

$$v(\mathbf{x},0)=u_t(\mathbf{x},0)=f(\mathbf{x},0)p(\mathbf{x}),$$

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

$$v_t(\mathbf{x},t)=\mathcal{A}v(\mathbf{x},t)+\frac{f_t(\mathbf{x},t)}{f(\mathbf{x},0)}v(\mathbf{x},0)$$

$$v(\mathbf{x},t)=0\text{mbox}and{\partial }_{\nu }v(\mathbf{x},t)=G_t(\mathbf{x},t)$$

$$Lv(\mathbf{x},t)=v_t(\mathbf{x},t)-\mathcal{A}v(\mathbf{x},t)-\frac{g_t(\mathbf{x},t)}{g(\mathbf{x},0)}v(\mathbf{x},0)$$

$$J_{ϵ}(v)={\int }_0^T{\int }_{\Omega }|Lv(\mathbf{x},t){|}^{2}d\mathbf{x}dt+ϵ\Vert v{\Vert }_{H^{2,1}(\Omega \times [0,1])}^2$$

$$H=\left\{v\in H^{2,1}(\Omega \times (0,T))\text{mbox}thatsatisfies\text{eqref}3.4\right\}$$

$$⟨L(w+\mathcal{E}),L\varphi {⟩}_{L^2(\Omega \times [0,T])}+ϵ⟨w,\varphi {⟩}_{H^{2,1}(\Omega \times [0,T])}=0,$$

$$\begin{gathered} ⟨Lw,L\varphi {⟩}_{L^2(\Omega \times [0,T])}+ϵ⟨w,\varphi {⟩}_{H^{2,1}(\Omega \times [0,T])} \\ =-⟨L\mathcal{E},L\varphi {⟩}_{L^2(\Omega \times [0,T])}-ϵ⟨\mathcal{E},\varphi {⟩}_{H^{2,1}(\Omega \times [0,T])}. \end{gathered}$$

$$v(\mathbf{x},t)=\sum _{n=1}^{\infty }{v}_n(\mathbf{x}){\Psi }_n(t)\text{mbox}forall(\mathbf{x},t)\in \Omega \times [0,T]$$

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

$$v(\mathbf{x},t)=\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(t)$$

$$v_t(\mathbf{x},t)=\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n^{\prime }(t)$$

$$\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n^{\prime }(t)=\sum _{n=1}^N\mathcal{A}{v}_n(\mathbf{x}){\Psi }_n(t)+\frac{f_t(\mathbf{x},t)}{f(\mathbf{x},0)}\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(0).$$

$$\begin{gathered} \sum _{n=1}^N{v}_n(\mathbf{x}){\int }_0^T{\Psi }_m(t){\Psi }_n^{\prime }(t)dt=\sum _{n=1}^N\mathcal{A}{v}_n(\mathbf{x}){\int }_0^T{\Psi }_m{\Psi }_n(t)dt \\ +\sum _{n=1}^N{v}_n(\mathbf{x}){\Psi }_n(0){\int }_0^T\frac{f_t(\mathbf{x},t)}{f(\mathbf{x},0)}{\Psi }_m(t)dt \end{gathered}$$

$$\mathcal{A}V=SV$$

$$S=\Big{(}\int_{0}^{T}(\Psi_{m}(t)\Psi^{\prime}_{n}(t)-\frac{f_{t}({\bf x},t)}{f({\bf x},0)}\Psi_{n}(0)\Psi_{m}(t)dt\Big{)}_{m,n=1}^{\infty}.$$

$$\Vert V{\Vert }_{H^1(\Omega {)}^N}\le C[\Vert V{\Vert }_{H^{1/2}(\partial \Omega {)}^N}+\Vert {\partial }_{\nu }V{\Vert }_{H^{-1/2}(\partial \Omega {)}^N}].$$

$$V(\mathbf{x})$$

$${\partial }_{\nu }V(\mathbf{x})$$

$$\displaystyle=\Big{(}G({\bf x},T)-G({\bf x},0)-\int_{0}^{T}G({\bf x},t)\Psi_{n}({\bf x},t)dt\Big{)}_{n=1}^{N}$$

$$\Vert V{\Vert }_{H^1(\Omega {)}^N}\le C[\Vert |G(\cdot ,T)|+|G(\cdot ,0)|{\Vert }_{H^{-1/2}(\partial \Omega )}+\Vert G{\Vert }_{H^{-1/2,1}(\partial \Omega \times [0,T])}].$$

$$\Vert v{\Vert }_{H^{2,1}(\Omega )}C[\Vert |G(\cdot ,T)|+|G(\cdot ,0)|{\Vert }_{H^{-1/2}(\partial \Omega )}+\Vert G{\Vert }_{H^{-1/2,1}(\partial \Omega \times [0,T])}].$$

$$\Vert p{\Vert }_{L^2(\Omega )}\le C[\Vert |G(\cdot ,T)|+|G(\cdot ,0)|{\Vert }_{H^{-1/2}(\partial \Omega )}+\Vert G{\Vert }_{H^{-1/2,1}(\partial \Omega \times [0,T])}].$$

$$\Omega =(-R,R{)}^2$$

$$\mathcal{G}=\left\{(x_i,y_j)=(-R+(i-1)d_{\mathbf{x}},-R+(j-1)d_{\mathbf{x}}):1\le i,j\le N+1\right\}$$

Full text

Key words. parabolic equation, inverse source problem, coefficient inverse problem, numerical method, quasi-reversibility method

AMS Classification 35R30, 35K20

1 Introduction

The area of inverse source problems has many applications and it, therefore, attracts the attention of the scientific community, see e.g., [14, 13, 15, 16, 17, 29, 32, 34, 35]. The solutions of inverse source problems can be used to directly detect the source even when the source is inactive after a certain time. Here, we name some examples. In the case of the parabolic equation, the problem plays an important role in identifying the pollution sources in a river or a lake [14]. In the case of elliptic equations, the inverse source problem has applications in electroencephalography [1, 13]. In the case that the data are generated by an acoustic source, the governing equation is the hyperbolic one and the problem addresses ultrasonics imaging and photoacoustic tomography [1, 13]. In this paper, we propose a numerical method to solve an inverse source problem for parabolic equations. This problem is the linearization of a nonlinear coefficient inverse problem. Therefore, we can use it to solve a coefficient inverse problem.

Let $\Omega$ be a bounded domain in $\mathbb{R}^{d}$, $d\geq 1$, with smooth boundary $\partial\Omega$. Let $c$ be a function in the class $C^{1}(\Omega)$. Consider the function $u=u({\bf x},t)$, ${\bf x}\in\Omega$, $t>0$ that is governed by the following initial value problem

[TABLE]

where $\mathcal{A}$ is an elliptic operator independent of the time and $f({\bf x},t)p({\bf x})$ is the source function. The aim of this paper is to solve the following inverse source problem.

Problem 1.1** (Inverse source problem for parabolic equations).**

Let $T$ be a positive number. Assume the function $f({\bf x},t)$, $({\bf x},t)\in\Omega\times[0,T]$, is known and $f({\bf x},0)$ does not vanish at any point ${\bf x}$ in $\Omega$. Determine the function $p({\bf x})$, ${\bf x}\in\Omega$, from the measurement of the following data

[TABLE]

for all ${\bf x}\in\partial\Omega$ and $t\in[0,T].$

The uniqueness of Problem 1.1 when the source function is a combination of some Dirac functions is confirmed in [14] and a numerical method to reconstruct this source is studied in [2]. We also draw the reader to the conditional stability in [20, 30]. In the case when the governing equation is the heat equation and the source function does not depend on the second variable, a reconstruction formula is provided in [32]. Another related problem is the inverse problem of reconstructing the initial condition for parabolic equation. This problem is very important and interesting, see [31, 36, 33, 27, 38] for theoretical results and numerical methods. In the current paper, we introduce the following approach to solve Problem 1.1. We derive from a governing equation a new equation involving only one unknown. The solution to that equation will directly provide the knowledge of the desired source function. However, that equation is not a standard partial differential equation. In fact, it involves the initial condition of itself. We prove the stability of the inverse source problem based on the projection of this equation on a finite dimensional space. A theory to solve this partial differential equation is not available yet. To solve this equation, we employ the quasi-reversibility method. This method was first introduced by Lattès and Lions [28]. 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., [4, 6, 7, 8, 11, 12, 18, 26, 21, 34]. We refer the reader to [22] for a survey on this method. The solutions of partial differential equations due to the quasi-reversibility method are called regularized solution in the theory of ill-posed problems [37]. The convergence of the regularized solution to the true one for three main types of partial differential equations is well-known [22]. Recently, in [34], the second author proved a Lipchitz convergence of quasi-reversibility method for the hyperbolic operator that involves Volterra integrals. The proof for a Lipchitz convergence of the quasi-reversibility method for the parabolic operator including the initial condition when this initial condition takes some particular forms will be proved in our near future publication.

An application of the inverse source problem in this paper is to solve a coefficient inverse problem for the heat equation. Given an initial guess of the coefficient, we show that our inverse source problem is a linear “perturbation” of that nonlinear coefficient inverse problem near that initial guess. Hence, by repeatedly solving our inverse source problem, we can obtain the solution to the coefficient inverse problem, see Section 6 for details. It is worth mentioning that the optimal control method to solve nonlinear coefficient inverse problem is widely used [5, 9, 10, 19, 39] which provide good numerical results with reasonable initial guesses. We also refer the reader to [3, 23] for the convexification method and numerical results in 1D.

The paper is organized as follows. We propose an algorithm to solve Problem 1.1 in Section 2. In section 3, we study the stability of Problem 1.1 in an approximation context. Next, in Section 4, we present the details about the implementation of our algorithm. In Section 5, we show some numerical solutions to the inverse source problem. In Section 6, we solve the nonlinear coefficient inverse problem from which the inverse source problem above arises. Section 7 is for concluding remarks.

2 The inversion method

Define the function

[TABLE]

Since $\mathcal{A}$ does not depend on $t$, it follows from the partial differential equation in (1.1) that

[TABLE]

for all ${\bf x}\in\Omega,$ $t\in(0,T).$ The initial condition for the function $v$ can be computed as

[TABLE]

which implies

[TABLE]

Substituting this into (2.2), we obtain

[TABLE]

for all ${\bf x}\in\Omega$, $t\in[0,T].$ Note that equation (2.4) does not depend on the function $p({\bf x})$.

Problem 1.1 becomes the problem of computing the function $v$ that satisfies (2.4) and the lateral Cauchy conditions

[TABLE]

for all ${\bf x}\in\partial\Omega,t\in[0,T].$

Remark 2.1**.**

We consider the function $G_{t}({\bf x},t)$ as our “indirect” data. In this paper, we test our method with noisy data $G_{t}({\bf x},t)=G_{t}({\bf x},t)(1+\delta(-1+2{\rm rand}))$ where $\delta$ is the noise level and ${\rm rand}$ is the uniformly distributed random number taking values in $[0,1].$ In this paper, $\delta=0\%,5\%$ and $10\%.$

Remark 2.2**.**

Our main idea when deriving (2.4) is that we want to eliminate one unknown so that we can arrive at the situation of one unknown and one equation. This strategy was applied in our research group in many publications; see e.g., [25, 35, 34]. Among them, the most similar idea to derive (2.4) is in [34] when the source term of a hyperbolic equation is eliminated. The main difference from (2.4) is that the corresponding equation in [34] is an integro-differential equation, which is not applicable in the current paper.

Assume that $v$ is known. Then, the desired function $p$ is computed via (2.3). However, due to the presence of the term $v({\bf x},0)$, equation (2.4), together with the lateral data in (2.5), is not a standard partial differential equation. A theortical method to solve it is not yet available. We solve (2.4) and (2.5) by the quasi-reversibility method. Define the operator

[TABLE]

for all function $v\in C^{2}(\overline{\Omega}\times[0,T])$. Given $\epsilon>0,$ we minimize the functional

[TABLE]

subject to the constraints in (2.5).

The following proposition guarantees that $J_{\epsilon}$ has a unique minimizer in $H$.

Proposition 2.1**.**

Assume that the set

[TABLE]

is nonempty. Then, for each $\epsilon>0,$ the function $J_{\epsilon}$ has a unique minimizer in $H$.

The proof of this proposition follows the proof of Proposition 3.1 in [35] for the time independent case. We present the proof for the time dependent case here for the connivence of the reader.

Proof of Proposition 2.1.

Let $\mathcal{E}$ be a function in $H$. Denote by $H_{0}$ the space $H-\mathcal{E}.$ Introduce $w=v-\mathcal{E}$. Then, minimizing $J_{\epsilon}(v)$ for $v$ in $H$ is equivalent to minimizing $J_{\epsilon}(w+\mathcal{E})$ for $w$ in $H_{0}.$ If $w\in H_{0}$ is a minimizer of $J_{\epsilon}(w+\mathcal{E})$ in $H_{0}$, then, by the variational principle,

[TABLE]

which is equivalent to

[TABLE]

The left hand side of (2.8) defines a new inner product $\{\cdot,\cdot\}$ in $H^{2,1}(\Omega\times[0,T])$. We have $\{w,w\}\geq\epsilon\|w\|_{H^{2,1}(\Omega\times[0,T])}^{2}$ and $\{w,w\}\leq C\|w\|_{H^{2,1}(\Omega\times[0,T])}^{2}$ for some constant $C$ due to the trace theory and the assumption that $\mathcal{A}$ is a second order elliptic operator. Hence, $\{\cdot,\cdot\}$ is equivalent to the standard inner product of $H^{2,1}(\Omega\times[0,T])$. On the other hand, the right hand side of (2.8) is a bounded linear operator defined on $H^{2,1}(\Omega\times[0,T])$. The existence and the uniqueness of a function $w$ satisfying (2.8) follows from the Riesz representation theorem. ∎

Remark 2.3**.**

*The unique minimizer of $J_{\epsilon}$ is call the regularized solution to (2.4) and (2.5). *

Our method to solve Problem 1.1 is summarized in Algorithm 1. In practice, we implement Algorithm 1 in the finite difference scheme. We present the implementation of Algorithm 1 with the finite difference method in the Section 4.

3 A Lipschitz estimate based on a truncation of the Fourier series

Let $\{\Psi_{n}\}_{n=1}^{\infty}$ be an orthonormal basis of $L^{2}(0,T).$ For each ${\bf x}\in\Omega$, we can write

[TABLE]

where $v({\bf x},t)$ is the function defined in (2.1). Here,

[TABLE]

Approximate the series in (3.1) by

[TABLE]

for $({\bf x},t)\in\Omega\times[0,T]$ for some number $N>0$. We also write

[TABLE]

Plugging (3.3) and (3.4) into (2.4), we have

[TABLE]

Multiply both side of the equation above by $\Psi_{m}(t)$ for each $m\in\{1,\dots,N\}$ and then integrate the resulting equation on $[0,T]$. We obtain

[TABLE]

for all ${\bf x}\in\Omega.$ Define $V({\bf x})=(v_{1}({\bf x}),\dots,v_{N}({\bf x}))$. It follows by (3.5) and the fact that $\Psi_{m}$ that the vector valued function $V({\bf x})$ satisfies the system

[TABLE]

where $S$ is a $d\times d$ matrix valued function given by

[TABLE]

Since $f$ is a smooth function, so is $S$. By a standard compact argument for elliptic equation, we can find a constant $C$ depending only on $\mathcal{A}$, $N$ and $\Omega$ such that

[TABLE]

It follows from (1.1), (1.2), (2.1) and (3.2) that

[TABLE]

on $\partial\Omega$. Hence, by (3.7),

[TABLE]

Using (3.3), we have

[TABLE]

As a result, using (2.3) and the trace theory, we get

[TABLE]

In summary, we have proved the following theorem.

Theorem 3.1**.**

Assume that the function $v({\bf x},t)=u_{t}({\bf x},t)$ is well-approximated by the Fourier sum as in (3.3) for some integer $N$ where $u({\bf x},t)$ is the solution to (1.1). Then, there exists a constant $C$ depending only on $\mathcal{A}$, $N$ and $\Omega$ such that

[TABLE]

Theorem 3.1 implies the Lipschitz stability for Problem 1.1 in the finite dimensional space spanned by $\{\Psi_{1},\dots,\Psi_{N}\}$. Studying the stability when $N$ tends to $\infty$ is extremely challenging and is out of the scope of this paper.

Remark 3.1**.**

The assumption about the well-approximation in Theorem 3.1 is verified numerically in some recent works of our research group. This verification for elliptic equation can be found in [35] and the one for parabolic equation is in [31]. In those papers, the basis $\{\Psi_{n}\}_{n=1}^{\infty}$ is taken from [24].

4 The finite difference method to find the regularized solution

In this section, the domain $\Omega$ is set to be a square in $\mathbb{R}^{2}$; i.e,

[TABLE]

where $R$ is a positive number. Let $N_{\bf x}$ and $N_{t}$ be positive integers. Set $d_{\bf x}=2R/N_{\bf x}$ and $d_{t}=T/Nt.$ We define a set of grid points on $\overline{\Omega}$

[TABLE]

and define a uniform partition on the time domain $[0,T]$ as

[TABLE]

For the simplicity in implementation, in this section, we modify the $H^{2,1}$ norm in the regularization term in (2.7) to the $H^{1}$ norm. In other words,

[TABLE]

for all $v\in H^{2,1}(\Omega\times[0,T]).$

Remark 4.1**.**

We replace the norm in regularization term $\epsilon\|v\|_{H^{2,1}(\Omega\times T)}^{2}$ by the $H^{1}$–norm because the $H^{1}$–norm is easier to implement. On the other hand, we have not observed any instabilities probably because the number $100\times 100$ of grid points we use is not too large and all norms in finite dimensional spaces are equivalent.

Remark 4.2** (The choice of $\epsilon$).**

We observe numerically that if $\epsilon$ is larger than $10^{-5}$, the reconstructed images of the source function are good but the reconstructed values are low and if $\epsilon<10^{-9}$, our method breaks down. We choose $\epsilon=10^{-8}$ in all our numerical tests. Note that this choice of $\epsilon$ is independent of the noise level, which is, in practice, supposed to be unknown.

The finite difference version of $J_{\epsilon},$ still named as $J_{\epsilon}$, reads

[TABLE]

Here, $L^{d_{\bf x},d_{t}}$ is the approximation of $L$ in the finite difference scheme and $\nabla^{d_{\bf x}}$ is the finite difference gradient. From now on, for the simplicity and to minimize the effort of writing computational code, we consider the case

[TABLE]

for some function $c$ in $L^{\infty}(\Omega).$ In this case,

[TABLE]

and

[TABLE]

for all $1\leq i,j\leq N_{{\bf x}}$ and $1\leq k\leq N_{t}+1.$ Introduce the $N$ dimensional vector $\mathfrak{v}$, $N=(N_{{\bf x}}+1)^{2}(N_{t}+1)$, whose $n^{\rm th}$ entry is given by

[TABLE]

where

[TABLE]

Then, we can rewrite (4.2) as

[TABLE]

where the $N\times N$ matrix $D$ is described as follows. For each $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ with $2\leq i,j\leq N_{{\bf x}}$ and $2\leq l\leq N_{t}$,

the $nn^{\rm th}$ entry $D_{nn}$ is given by $\frac{1}{d_{t}}+\frac{4}{d_{{\bf x}}^{2}}-c(x_{i},y_{j})$; 2. 2.

the $nm^{\rm th}$ entry $D_{nm}$ is given by $-\frac{1}{d_{t}}$ where $m=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l-1$ for $3\leq l\leq N_{t}$; 3. 3.

the $nm^{\rm th}$ entry $D_{nm}$ is given by $-\frac{1}{d_{t}}-\frac{f_{t}(x_{i},y_{j},t_{l})}{f(x_{i},y_{j},t_{1})}$ where $m=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l-1$ for $l=2$; 4. 4.

the $nm^{\rm th}$ entry $D_{nm}$ is given by $-\frac{1}{d_{{\bf x}}^{2}}$ where $m=(i\pm 1-1)(N_{{\bf x}}+1)(N_{t}+1)+(j\pm 1-1)(N_{t}+1)+l-1$ for $2\leq l\leq N_{t}$; 5. 5.

the other entries are [math].

We next define the matrices $D_{x}$ and $D_{y}$ such that $(D_{x}\mathfrak{v},D_{y}\mathfrak{v})=\nabla^{d_{\bf x}}v$. For each $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ with $2\leq i,j\leq N_{{\bf x}}+1$ and $1\leq l\leq N_{t}+1$,

the $nn^{\rm th}$ entry of $D_{x}$ and $D_{y}$ is given by $\frac{1}{d_{\bf x}}$; 2. 2.

the $nm^{\rm th}$ entry of $D_{x}$ is given by $-\frac{1}{d_{\bf x}}$ for $m=(i-1-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$; 3. 3.

the $nm^{\rm th}$ entry of $D_{y}$ is given by $-\frac{1}{d_{\bf x}}$ for $m=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1-1)(N_{t}+1)+l$; 4. 4.

other entries are $0.$

The finite difference version of $J_{\epsilon}$, defined in (4.1), becomes

[TABLE]

Hence, due to (4.3), since $v$ is a minimizer of $J_{\epsilon}$, $\mathfrak{v}$ satisfies the equation

[TABLE]

We next consider the boundary conditions for $v$ in (2.5). In the finite difference scheme, the first condition in (2.5) reads for $l=1,2,\dots,N_{t}+1,$

[TABLE]

for all $i\in\{1,N_{{\bf x}}+1\}$ and $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ or $i\in\{1,\dots,N_{{\bf x}}+1\}$ and $j\in\{1,N_{{\bf x}}+1\}$. Therefore, due to (4.3), we can write this condition as

[TABLE]

where $K_{1}$ is defined as follows. For $l\in\{1,2,\dots,N_{t}+1\},$

the $nn^{\rm th}$ entry of $K_{1}$ is set to be $1$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for some $i\in\{1,N_{\bf x}+1\}$, $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ or $i\in\{1,2,\dots,N_{\bf x}+1\}$, $j\in\{1,N_{{\bf x}}+1\}$; 2. 2.

the other entries of $K_{1}$ are [math].

The second condition in (2.5) is rewritten as

[TABLE]

where the vector $\mathfrak{g}$ is the lineup version of the data $G_{t}$

[TABLE]

for all $i\in\{1,N_{{\bf x}}+1\}$ and $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ or $i\in\{1,\dots,N_{{\bf x}}+1\}$ and $j\in\{1,N_{{\bf x}}+1\}$ and the matrix $K_{2}$ is defined as follows. For all $l\in\{1,2,\dots,N_{t}+1\}$,

the $nn^{\rm th}$ entry of $K_{2}$ is $\frac{1}{d_{\bf x}}$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for $i\in\{1,N_{\bf x}+1\}$, $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ or $i\in\{1,2,\dots,N_{\bf x}+1\}$, $j\in\{1,N_{{\bf x}}+1\}$; 2. 2.

the $nm^{\rm th}$ entry of $K_{2}$ is $-\frac{1}{d_{\bf x}}$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for $i=1$, $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ and $m=(i+1-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$; 3. 3.

the $nm^{\rm th}$ entry of $K_{2}$ is $-\frac{1}{d_{\bf x}}$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for $i=N_{\bf x}+1$, $j\in\{1,2,\dots,N_{{\bf x}}+1\}$ and $m=(i-1-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$; 4. 4.

the $nm^{\rm th}$ entry of $K_{2}$ is $-\frac{1}{d_{\bf x}}$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for $i\in\{2,\dots,N_{{\bf x}}\}$, $j=1$ and $m=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j+1-1)(N_{t}+1)+l$; 5. 5.

the $nm^{\rm th}$ entry of $K_{2}$ is $-\frac{1}{d_{\bf x}}$ if $n=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1)(N_{t}+1)+l$ for $i\in\{2,\dots,N_{{\bf x}}\}$, $j=N_{\bf x}+1$ and $m=(i-1)(N_{{\bf x}}+1)(N_{t}+1)+(j-1-1)(N_{t}+1)+l$; 6. 6.

the other entries of $K_{2}$ are 0.

Combining (4.5), (4.6) and (4.7), we obtain

[TABLE]

Since $\epsilon$ is a small number, it is acceptable that we modify the equation above by a more “stable” one

[TABLE]

The analysis in this section is summarized in the following proposition.

Proposition 4.1**.**

The source function $p({\bf x})$ can be computed by

solve (4.8) for $\mathfrak{v}$, the “line up” version of $v$; 2. 2.

compute the function $v_{\rm comp}$ using (4.3); 3. 3.

calculate $p_{\rm comp}({\bf x})=\frac{v_{\rm comp}({\bf x})}{f({\bf x},0)}$.

5 Numerical results

We test our numerical method when $R=1$ and $\Omega,$ therefore, is $(-1,1)^{2}.$ Also, we choose $T=0.2$, see the Remark 1 for this choice of $T$.

Remark 5.1** (Choose $T$).**

We numerically choose $T$ by examining the $L^{2}$ norm of the data $G_{t}({\bf x},t)$ on $\partial\Omega$ as a function in $T$. Define

[TABLE]

The graph of the function $\gamma$ is displayed in Figure 1, showing that the data is largest on $[0,0.2]$. This means the data contains most important information about the source in this interval. We therefore choose $T=0.2$ for all of our numerical tests.

We chose $N_{{\bf x}}=100$ and $N_{t}=60$ in this section. In all tests, the known function $f$ is chosen as

[TABLE]

and the known function $c({\bf x})$ is set to be

[TABLE]

In this section, we show the following numerical results.

Test 1. In this test, the true source function $p_{\rm true}$ is smooth and given by

[TABLE]

The numerical result for this test is displayed in Figure 2.

It is evident that our method well reconstructs the source function $p_{\rm true}$. The location and shape of the circular “inclusion” can be identified. The true maximum value of the inclusion is $1$. The reconstructed maximum value of the inclusion is computed with small errors. In fact,

when $\delta=0\%,$ $\max_{{\bf x}\in\Omega}p_{\rm comp}({\bf x})=0.991$ and the coresponding relative error is $0.9\%$; 2. 2.

when $\delta=5\%,$ $\max_{{\bf x}\in\Omega}p_{\rm comp}({\bf x})=0.976$ and the coresponding relative error is $2.4\%$; 3. 3.

when $\delta=10\%,$ $\max_{{\bf x}\in\Omega}p_{\rm comp}({\bf x})=1.048$ and the coresponding relative error is $4.8\%$.

Test 2. We test our method for the case when $p_{\rm true}$ is given by the smooth function

[TABLE]

In this test, the true source function has a negative “inclusion” and a positive one. The numerical results for this test are displayed in Figure 3.

The true and computed local extreme values of the source function at two inclusions are displayed in Table 1. This table show that our method is stable with respect to noise.

Test 3. We next check the case when the source function is not smooth. In this case, we consider the piecewise constant function

[TABLE]

The graph of the function $p_{\rm true}$ has two “inclusions” with different shapes, a disk and an ellipse. The graphs of the true and computed source function are displayed in Figure 4.

The reconstruction of the image of the source function in this test is acceptable. Table 2 shows the strength of our method in the sense that we can reconstruct the values of those two inclusions with acceptable error.

Test 4. We test the nonsmooth true source function again with a complicated support. The true function is the characteristic function of the letter $\Omega.$ The numerical results for this test are shown in Figure 5.

Note that the maximal value of the reconstructed functions are acceptable. When $\delta=0\%,$ $\max p_{\rm comp}=0.97$ (relative error 3%). When $\delta=5\%,$ $\max p_{\rm comp}=0.97$ (relative error 3%). When $\delta=10\%,$ $\max p_{\rm comp}=1.073$ (relative error 7.3%).

Remark 5.2**.**

Despite of the presence of the initial condition in equation (2.4), it is evident that the quasi-reversibility method provides good numerical results with small relative errors.

6 Application to a coefficient inverse problem

In this section, we propose a numerical method to solve a severely ill-posed nonlinear coefficient inverse problem.

6.1 The problem statement

Problem 1.1 arises from a coefficient inverse problem for parabolic equations. For the simplicity, consider the problem of determining the coefficient $c({\bf x})$ from the measurements of $\partial_{n}\mathfrak{u}({\bf x},t)$ on $\partial\Omega\times[0,T]$ where, $\mathfrak{u}({\bf x},t)$ is the solution of the following problem

[TABLE]

Assume that the initial condition $g({\bf x})>0$ for all ${\bf x}\in\Omega$ and the boundary condition $g_{1}$ satisfying $g_{1}({\bf x},0)=g({\bf x})$ for all ${\bf x}$ in $\partial\Omega.$ Consider the following nonlinear inverse problem.

Problem 6.1**.**

Let $T>0.$ Determine the coefficient $c({\bf x})$, ${\bf x}\in\Omega$ from the measurement of

[TABLE]

for all ${\bf x}\in\partial\Omega,$ $t\in[0,T].$

Problem 6.1 and its related versions are 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., [5, 9, 10, 19, 39] and references therein. The main drawback of this method is that the initial guess for the true solution is important to obtain numerical results. Un like this, we assume that we do not have any advanced knowledge of the true solution to Problem 6.1 and take the initial guess as a constant function.

Consider the circumstance that an initial guess for the function $c$, named as $c_{0}$, is known. Then, we write

[TABLE]

Denote by the function $\mathfrak{u}_{0}({\bf x})$ the solution of (6.1) with $c_{0}$ replacing $c$ and let $w=\mathfrak{u}-\mathfrak{u}_{0}$. It is not hard to see that

[TABLE]

Since $c_{0}$ is an initial guess of $c$, we can replace the function $\mathfrak{u}$ in the differential equation in (6.3) above by $\mathfrak{u}_{0}$ to obtain

[TABLE]

which leads to a particular case of Problem 1.1 with $f=\mathfrak{u}_{0}$. We can compute $p({\bf x})$ and therefore $c({\bf x})$ via solving Problem 1.1 for the heat equation (6.4). Denoting the computed $c({\bf x})$ by $c_{1}({\bf x})$ and let $\mathfrak{u}_{1}({\bf x},t)$ be the solution to (6.1) with $c=c_{1}$. We then find $c_{2}$ by solving Problem 1.1 for the heat equation (6.4) with $\mathfrak{u}_{1}$ replacing $\mathfrak{u}_{0}$. The process is repeated to compute $c_{3},c_{4},\dots$ and we choose $c_{\rm comp}=c_{n^{*}}$ when $n^{*}$ is a fixed positive integer. We summarize this numerical method to compute $c$ in Algorithm 2.

Remark 6.1**.**

Imposing assumption (6.2), where the function $c_{0}$ is known and the unknown function $p$ is small, only plays the role of the suggestion for the “linearization” analysis. However, in the reverse direction, the numerical results show that Algorithm 2 can be applied and provide good numerical results even in the case when $c_{0}$ is far away from the function $c$. Here, we understand by “$c_{0}$ is far away from the function $c$” in two senses: (1) the complicated geometry of the true function $c$ and (2) the high contrasts.

Remark 6.2**.**

The difficulty about the presence of an initial guess $c_{0}$ can be overcome using the quasi-reversibility method. The authors of [3, 23] introduce a convex functional, which minimizer yields the solution of the problem under consideration, by combining the quasi-reversibility method and the Carleman weight functions. Numerical results in 1D are presented in [3]. It is important and interested to numerically test their method in higher dimensions.

In the next subsection, we will show some numerical results. We also display the graph of the relative difference

[TABLE]

to show the convergence of Algorithm 2.

6.2 Numerical results

We perform two numerical results due to Algorithm 2 below. In these tests, the noise level is $5\%$. The background function $c_{0}({\bf x})=1$ for all ${\bf x}\in\Omega.$

Test 5. The function $c_{\rm true}$ is the step function taking value 3 inside a letter $\Sigma$ and $1$ otherwise. We display the obtained numerical results in Figure 6.

The reconstructed image “$\Sigma$” meets the expectation although the initial guess $c_{0}=1$ is far away from the true function $c_{\rm true}$. The true maximal value of the function $c_{\rm true}$ is 3 and the reconstructed one is 3.56. The relative error is $18\%.$

Test 6. In test 6, the function $c_{\rm true}$ is given by

[TABLE]

The image of the function $c_{\rm true}$ has a horizontal rectangle and a vertical rectangle. Due to the geometry and the high value, $c_{\rm true}$ is far away from the background $c_{0}=1$. We display the obtained numerical results in Figure 7. Despite of the “bad” initial guess $c_{0}$, the rectangles can be seen after a few iterations. We note that the reconstructed images and value are improved with more iterations. The true maximum value of $c_{\rm true}$ is 5 and the reconstructed one is 5.94. The relative error is 18.8%.

Remark 6.3**.**

We observe that with the choice of $c_{0}$ as the background constant, the first reconstructed function $c_{1}$ is poor. Then, in the next two iterations, the quality of the reconstructed function improves significantly. Figures 6f and 7f show that the sequence $\{c_{n}\}_{n\geq 1}$ converges at the very fast rate.

7 Concluding remarks

In this paper, we have proposed a method to solve an inverse source problem for parabolic equations. The stability of this problem is proved in an approximation context. To compute the numerical solutions to this inverse source problem, we derived an equation whose solution directly provides the desired solution of our inverse source problem. However, this equation is not a standard parabolic equation. A theory to solve it is not yet available. We therefore employ the quasi-reversibility method to find its solution. Since the inverse source problem in this paper is a linearization of a nonlinear coefficient inverse problem, we use the proposed method to establish an iterative method to solve that nonlinear coefficient inverse problem. Numerical results were presented.

Acknowledgments

The authors sincerely appreciate Michael V. Klibanov for many fruitful discussions. The work of the second author was supported by US Army Research Laboratory and US Army Research Office grant W911NF-19-1-0044.