Numerical Solutions to the Robin Inverse Problem with Nonnegativity Constraints

Weifu Fang¹, Fu-Rong Lin²

¹Department of Mathematics and Statistics, Wright State University, Dayton, USA.
²Department of Mathematics, Shantou University, Shantou, China.

DOI: 10.4236/jamp.2022.106137                

Abstract

We present iterative numerical methods for solving the inverse problem of recovering the nonnegative Robin coefficient from partial boundary measurement of the solution to the Laplace equation. Based on the boundary integral equation formulation of the problem, nonnegativity constraints in the form of a penalty term are incorporated conveniently into least-squares iteration schemes for solving the inverse problem. Numerical implementation and examples are presented to illustrate the effectiveness of this strategy in improving recovery results.

Keywords

Robin Inverse Problem, Ill-Posed Problem, Laplace Equation, Boundary Inte-gral Equation, Nonnegativity Constraint, Penalty Method

1. Introduction

Consider the Robin boundary value problem for the Laplace equation on a smooth bounded domain $\Omega \subset R^2$ :

$\{\begin{array}{ll}\Delta U=0\hfill & \text{in}\Omega ,\hfill \\ \frac{\partial U}{\partial \nu }+pU=g\hfill & \text{on}\partial \Omega =\Gamma .\hfill \end{array}$ (1.1)

Here $p=p\left(x\right)\ge 0$ is the Robin coefficient and $g=g\left(x\right)\ge 0$ is a prescribed input function on $\Gamma$. It is assumed that the support of p (denoted by ${\Gamma }_1$ ) and the support of g are nonempty but non-overlapping. Then the Robin inverse problem is described as: to recover the Robin coefficient p on ${\Gamma }_1$ from given values of the solution $U=u_0$ on ${\Gamma }_0\subset \Gamma$ with ${\Gamma }_0\cap {\Gamma }_1=\varnothing$. Such problems arise from the modeling of various nondestructive evaluation techniques where the solution U can be measured directly on an accessible part ${\Gamma }_0$ of the boundary in hope to extract the quantity of interest p defined on an inaccessible portion ${\Gamma }_1$ of the boundary. For example, in the evaluation of metal-to-silicon contact quality in semiconductor transistors (e.g. [1] [2] [3]), voltage measurements of U on an accessible ${\Gamma }_0$ in response to an applied current input pattern g are used to infer the contact quality modeled by p on the inaccessible contact window ${\Gamma }_1$. A similar situation arises from detection of material corrosion damages (e.g. [4] [5]), where static potential measurement u on an accessible part of the material boundary ${\Gamma }_0$ is taken to detect possible corrosion characterized by the conductivity profile p on an inaccessible $\Gamma$.

For this inverse problem, there have been many numerical studies in the literature, based either on the PDE model (e.g. [5] [6] [7] and references therein) or on an integral equation formulation (e.g. [3] [8] [9] and references therein). Most of these papers have focused on the recovery of the Robin coefficients by iterative schemes for the minimizers of some objective functional that consists of both a data fitting term and a certain form of regularization, and the nonnegativity constraint for the Robin coefficient is often ignored, except for some simple truncation strategies on the side to ensure the well-posedness of the forward solver. It is well known that the inclusion of nonnegativity constraints improves the quality of solutions to ill-posed inverse problems when the quantity of interest is known to be nonnegative (e.g. [10], Chapter 9). In this study, we propose a strategy to incorporate nonnegativity constraints in the iterative schemes, in the form of an additional term in the objective functional to penalize negativity throughout the iteration process. This approach is natural in the least-squares formulation for solving inverse problems or illposed problems in general and provides a convenient and computationally economical alternative to any constrained optimization methods. We demonstrate the simplicity of this strategy in both linear and nonlinear least-squares solution methods for the Robin inverse problem, and through numerical examples, we illustrate the effectiveness of the proposed approach in dealing with the ill-posedness, resulting in much improved recovery results.

Our presentation is organized as follows. In Section 2, we reformulate the problem (1.1) into its equivalent boundary integral equation formulation. Based on this formulation, we present two least-squares solution methods for the inverse problem in Section 3, with nonnegativity constraints incorporated into the iterative schemes. In Section 4, we conclude with discussions of numerical implementation of the recovery algorithms, and examples of numerical results to illustrate how the added nonnegativity constraints have helped to improve solutions to the inverse problem.

2. Integral Equation Formulation

Because the equation in (1.1) is Laplacian over the domain $\Omega$ while the quantities of interest, the Robin coefficient p and data measurement $u_0$, are both defined on the boundary $\Gamma$ only, it is natural and advantageous to formulate this boundary value problem (1.1) as a boundary integral equation on $\Gamma$.

Let $\Phi =\Phi \left(x\mathrm{,}y\right)$ stand for the fundamental solution for the Laplacian in $R^2$ :

$\Phi \left(x,y\right)=\frac{1}{2\pi }\ln \frac{1}{\left|x-y\right|}\text{for}x\ne y.$

Then the boundary value problem (1.1) for $U\in H^1\left(\Omega \right)$ in $\Omega$ is equivalent to the following boundary integral equation for $u\in H^{1/2}\left(\Gamma \right)$ on $\Gamma$ (as the trace of U) (see e.g. [11] [12] [13]):

$\frac{1}{2}u\left(x\right)+{\displaystyle {\int }_{\Gamma }}\left(\frac{\partial \Phi \left(x\mathrm{,}y\right)}{\partial {\nu }_y}+p\left(y\right)\Phi \left(x\mathrm{,}y\right)\right)u\left(y\right)\text{d}{s}_y={\displaystyle {\int }_{\Gamma }}\Phi \left(x\mathrm{,}y\right)g\left(y\right)\text{d}{s}_y\mathrm{,}x\in \Gamma \mathrm{.}$ (2.1)

In the operator form, (2.1) can be written as

$\left(\frac{1}{2}\mathcal{I}+\mathcal{D}\right)u+\mathcal{S}\left(pu\right)=\mathcal{S}g\mathrm{,}$ (2.2)

with the single and double-layer potential operators defined by

$\left(\mathcal{S}u\right)\left(x\right)={\displaystyle {\int }_{\Gamma }}\Phi \left(x,y\right)u\left(y\right)\text{d}{s}_y\text{and}\left(\mathcal{D}u\right)\left(x\right)={\displaystyle {\int }_{\Gamma }}\frac{\partial \Phi \left(x,y\right)}{\partial {\nu }_y}u\left(y\right)\text{d}{s}_y\text{for}x\in \Gamma .$

Note that the operators have the following mapping properties (e.g. [13]): $\mathcal{S}:H^{-1/2}\left(\Gamma \right)\to H^{1/2}\left(\Gamma \right)$ and $\mathcal{D}:H^{1/2}\left(\Gamma \right)\to H^{1/2}\left(\Gamma \right)$. It is well known that the Robin boundary-value problem (2.1) admits unique solution $U\ge 0$ in $\Omega$ for given nonnegative model parameters p and g on $\Gamma$, and thus the equivalent boundary integral equation (2.1) yields nonnegative solution u on $\Gamma$ as well (see e.g. [11]).

Alternatively, we can use single-layer potential to express the solution $U\left(x\right)$ on $\Omega$ by

$U\left(x\right)={\displaystyle {\int }_{\Gamma }}\Phi \left(x\mathrm{,}y\right)\phi \left(y\right)\text{d}{s}_y\mathrm{,}x\in \Omega \mathrm{,}$ (2.3)

where the density function $\phi$ on $\Gamma$ solves the boundary integral equation

$\frac{1}{2}\phi \left(x\right)+{\displaystyle {\int }_{\Gamma }}\frac{\partial \Phi \left(x\mathrm{,}y\right)}{\partial {\nu }_x}\phi \left(y\right)\text{d}{s}_y+p\left(x\right){\displaystyle {\int }_{\Gamma }}\Phi \left(x\mathrm{,}y\right)\phi \left(y\right)\text{d}{s}_y=g\left(x\right)\mathrm{,}x\in \Gamma \mathrm{.}$ (2.4)

In operator form, (2.4) becomes

$\left(\frac{1}{2}\mathcal{I}+{\mathcal{D}}^{\prime }\right)\phi +p\cdot \mathcal{S}\phi =g,$

where the dual operator ${\mathcal{D}}^{\prime }$ of $\mathcal{D}$ is given by

$\left({\mathcal{D}}^{\prime }\phi \right)\left(x\right)={\displaystyle {\int }_{\Gamma }}\frac{\partial \Phi \left(x\mathrm{,}y\right)}{\partial {\nu }_x}\phi \left(y\right)\text{d}{s}_y\mathrm{,}x\in \Gamma \mathrm{.}$

In our numerical implementation, we use the direct formulation (2.1) as the model equation that relates between u and p, and only use the alternative formulation (2.3) - (2.4) for generating synthetic data $u_0$ which, after the addition of random noise, are used to test our numerical algorithms for recovering the Robin coefficient p.

3. Solution Methods for the Inverse Problem

A common approach to the inverse problem of finding p on ${\Gamma }_1$ from given boundary measurement $u_0$ of u on ${\Gamma }_0$ is to regard u as dependent on p through the model equation (2.2), and then solve the data-fitting nonlinear equation $u=u_0$ for p by least-squares methods (see e.g. [3] [7]). In this section, we consider two slightly different solution methods, each of which is based on solving $\left(u\mathrm{,}p\right)$ simultaneously from the system consisting of both the model equation (2.2) and the data fitting equation $u=u_0$ on ${\Gamma }_0$ by least-squares (see e.g. [9] [11] [14]). Then in the iterative scheme for solving the least-squares problem, we will introduce nonegativity constraints in the form of a penalty term in the minimizing objective functional, so to steer toward nonnegative solutions.

3.1. A Linear Least-Squares Method

As in [9] [11], we introduce a new variable

$v\left(x\right)=p\left(x\right)u\left(x\right)\mathrm{,}x\in {\Gamma }_1$ (3.1)

in place of $p\left(x\right)$. Then equation (2.2) becomes linear in both $u$ and $v$ :

$\left(\frac{1}{2}\mathcal{I}+\mathcal{D}\right)u+{\mathcal{S}}_{1}v=\mathcal{S}g,$ (3.2)

where

$\left({\mathcal{S}}_{1}v\right)\left(x\right)={\displaystyle {\int }_{{\Gamma }_1}}\Phi \left(x\mathrm{,}y\right)v\left(y\right)\text{d}{s}_y\text{for}x\in \Gamma \mathrm{.}$

Denote the restriction operator from $\Gamma$ to ${\Gamma }_0$ by ${\mathcal{R}}_0\mathrm{<mo>\; :\; </mo>}{L}^2\left(\Gamma \right)\to L^2\left({\Gamma }_0\right)$. That is, $\left({\mathcal{R}}_{0}u\right)\left(x\right)=u\left(x\right)$ for $x\in {\Gamma }_0$. Then the given measurement $u_0$ of u on ${\Gamma }_0$ can be expressed as

${\mathcal{R}}_{0}u=u_0.$ (3.3)

We recast the inverse problem of finding p from given $u_0$ as directly solving for $w={\left(u\mathrm{,}v\right)}^{\text{T}}\in L^2\left(\Gamma \right)\times L^2\left({\Gamma }_1\right)$ from (3.2) - (3.3) as a linear system:

$\left[\begin{array}{cc}\frac{1}{2}\mathcal{I}+\mathcal{D}& {\mathcal{S}}_1\\ {\mathcal{R}}_0& \mathcal{O}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}\mathcal{S}g\\ u_0\end{array}\right]\text{or}\mathcal{A}w=f\mathrm{.}$ (3.4)

Here $\mathcal{O}$ denotes the zero operator from $L^2\left({\Gamma }_1\right)$ to $L^2\left({\Gamma }_0\right)$. Once u on $\Gamma$ and v on ${\Gamma }_1$ are found from (3.4), then the Robin coefficient p on ${\Gamma }_1$ can be found from u and v by the simple relation (3.1).

System (3.4) is a linear system for $w={\left(u\mathrm{,}v\right)}^{\text{T}}$, but is ill-posed. In addition to using the classical Tikhonov regularization method to address the ill-posedness (as in [11]), we wish to incorporate an additional term to penalize negativity in w and encourage nonnegagive solutions. Thus we seek the minimizer of the following objective functional that consists of the system residual from (3.4), a regularization for w, and a penalty for negative part of w:

$J\left(w\right)=\frac{1}{2}{\Vert \mathcal{A}w-f\Vert }_2^2+\frac{\alpha }{2}\langle \mathcal{H}w\mathrm{,}w\rangle +\frac{\beta }{2}{\Vert w^{-}\Vert }_2^2\mathrm{,}$ (3.5)

with regularization parameter $\alpha >0$ and penalty parameter $\beta >0$. The regularization is defined through the operator

$\mathcal{H}w={\left(-D_p^{2}u\mathrm{,}-D_0^{2}v\right)}^{\text{T}}$

where $D_p^2$ and $D_0^2$ are the second-order derivative operator for periodic boundary condition and for zero boundary condition respectively, and the penalty for the negative part of w by

$w^{-}=-\min \left(0,w\right)=\frac{1}{2}\left(\left|w\right|-w\right)\equiv \psi \left(w\right).$

The function $\psi \left(s\right)$ is differentiable except at $s=0$, and when its derivative is needed, we use the derivative of its smooth approximation:

${{\psi }^{\prime }}_{\epsilon }\left(s\right)=\frac{1}{2}\left(\frac{s}{\sqrt{s^2+\epsilon }}-1\right)\text{where}{\psi }_{\epsilon }\left(s\right)=\frac{1}{2}\left(\sqrt{s^2+\epsilon }-s\right)$

for some very small $\epsilon >0$ (e.g. the machine epsilon). Hence we can arrive at the first-order optimality condition of $J\left(w\right)$ :

$\left({\mathcal{A}}^{\prime }\mathcal{A}+\alpha \mathcal{H}\right)w+\beta \psi \left(w\right){{\psi }^{\prime }}_{\epsilon }\left(w\right)={\mathcal{A}}^{\prime }f$ (3.6)

where ${\mathcal{A}}^{\prime }$ denotes the dual of $\mathcal{A}$.

In (3.6), while the parts from the linear system and Tikhonov regularization are linear, the part from the penalty term is not. Hence we devise an iterative scheme to solve (3.6) using the Gauss-Newton direction as follows. At each iterate $w^{\left(k\right)}$, the step increment $d^{\left(k\right)}$ is determined by the linear equation

$\left({\mathcal{A}}^{\prime }\mathcal{A}+\alpha \mathcal{H}+\beta {\left({{\psi }^{\prime }}_{\epsilon }\left(w^{\left(k\right)}\right)\right)}^2\right)d^{\left(k\right)}=-r^{\left(k\right)}$ (3.7)

where

$r^{\left(k\right)}={\mathcal{A}}^{\prime }\left(\mathcal{A}{w}^{\left(k\right)}-f\right)+\alpha \mathcal{H}{w}^{\left(k\right)}+\beta \psi \left(w^{\left(k\right)}\right){{\psi }^{\prime }}_{\epsilon }\left(w^{\left(k\right)}\right),$

and the next iterate is set to be

$w^{\left(k+1\right)}=w^{\left(k\right)}+d^{\left(k\right)}.$ (3.8)

We note that, when $\beta =0$, i.e. when the penalty term is not present, no iteration is needed, and this reduces to the direct method introduced in [11]. However, such a direct method may produce w that admits negative values, and thus it requires some remedy, such as simple truncation, in order to obtain the nonnegative Robin coefficient p from (3.1) (see [9] [11]). By incorporating the penalty term in the least-squares formulation, we are able to ensure the nonnegativity of the solution it produces, with only a few additional iterations, so the Robin coefficient p can be readily obtained from (3.1). Moreover, when the penalty term is active, it steers the iteration to produce better solutions, especially in situations where the ill-posedness is more severe.

3.2. A Nonlinear Least-Squares Method

In this method, we regard both u on $\Gamma$ and p on its support ${\Gamma }_1$ as variables to be found from given data $u_0$, and cast the inverse Robin problem as finding $z={\left(u,p\right)}^{\text{T}}\in L^2\left(\Gamma \right)\times L^2\left({\Gamma }_1\right)$ from the nonlinear system of equations consisting of both the model equation (2.2) and data fitting equation (3.3):

$\{\begin{array}{l}\left(\frac{1}{2}\mathcal{I}+\mathcal{D}\right)u+\mathcal{S}\left(pu\right)=\mathcal{S}g,\\ {\mathcal{R}}_{0}u=u_0,\end{array}\text{or}\mathcal{B}\left(z\right)=f.$ (3.9)

Here p in $z={\left(u\mathrm{,}p\right)}^{\text{T}}$ to be solved refers to the part of the Robin coefficient on its support ${\Gamma }_1$, while in the formulations p remains as defined on $\Gamma$ with the understanding that its support is contained in ${\Gamma }_1$. It is noted that $\mathcal{B}$ is linear in u and p individually, but not jointly in z.

To solve the nonlinear system (3.9) for $z={\left(u\mathrm{,}p\right)}^{\text{T}}$, we consider minimization of the least-squares functional of the system residual, with the addition of both a regularization term and a penalty term:

$J\left(z\right)=\frac{1}{2}{\Vert \mathcal{B}\left(z\right)-f\Vert }_2^2+\frac{\alpha }{2}\langle \mathcal{H}z\mathrm{,}z\rangle +\frac{\beta }{2}{\Vert z^{-}\Vert }_2^2$ (3.10)

with regularization $\mathcal{H}$ and penalty $z^{-}=\psi \left(z\right)$ as in (3.5). Note that the case without penalty (when $\beta =0$ ) was considered in [14]. For minimization of J above, we propose the Gauss-Newton iteration scheme: at each iterate $z^{\left(k\right)}$, the new increment $d^{\left(k\right)}$ is determined by the Gauss-Newton step from J:

$\left(\stackrel{\dot{}}{\mathcal{B}}{\left(z^{\left(k\right)}\right)}^{\prime }\stackrel{\dot{}}{\mathcal{B}}\left(z^{\left(k\right)}\right)+\alpha \mathcal{H}+\beta {\left({{\psi }^{\prime }}_{\epsilon }\left(z^{\left(k\right)}\right)\right)}^2\right)d^{\left(k\right)}=-r^{\left(k\right)}$ (3.11)

where

$r^{\left(k\right)}=\stackrel{\dot{}}{\mathcal{B}}{\left(z^{\left(k\right)}\right)}^{\prime }\left(\stackrel{\dot{}}{\mathcal{B}}\left(z^{\left(k\right)}\right)-f\right)+\alpha \mathcal{H}{z}^{\left(k\right)}+\beta \psi \left(z^{\left(k\right)}\right){{\psi }^{\prime }}_{\epsilon }\left(z^{\left(k\right)}\right),$

and $\stackrel{\dot{}}{\mathcal{B}}\left(z\right)$ is the derivative operator of $\mathcal{B}\left(z\right)$ given by

$\stackrel{\dot{}}{\mathcal{B}}\left(z\right)=\left[\begin{array}{cc}\frac{1}{2}\mathcal{I}+\mathcal{D}+\mathcal{S}\circ \left(p\mathcal{I}\right)& \mathcal{S}\circ \left(u\mathcal{I}\right)\\ {\mathcal{R}}_0& \mathcal{O}\end{array}\right]$

where $\circ$ denotes operator composition, and $\stackrel{\dot{}}{\mathcal{B}}{\left(z\right)}^{\prime }$ the dual of $\stackrel{\dot{}}{\mathcal{B}}\left(z\right)$. Once $d^{\left(k\right)}$ is found from (3.11), the next iterate is set to be

$z^{\left(k+1\right)}=z^{\left(k\right)}+d^{\left(k\right)}\mathrm{.}$ (3.12)

Note that the addition of the penalty term in (3.10) results in a similar iterative scheme (3.11), and it does not add much computational cost. As it steers the iterations to satisfy the nonnegativity constraints, it does improve the quality of the solutions, as numerical examples in the next section will illustrate.

4. Numerical Examples

Given a regular 1-periodic parametrization with counterclockwise orientation for the boundary $\Gamma$ :

$x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right)\right),0\le t\le 1,$

where $x_1\left(t\right),x_2\left(t\right)\in C_p^2\left[\mathrm{0,1}\right]$ and $\left|x^{\prime }\left(t\right)\right|>0$ for $0\le t\le 1$, the integral operators in (2.1) and (2.4) can be expressed explicitly as integral operators on $t\in \left[\mathrm{0,1}\right]$ for $u\left(t\right)=u\left(x\left(t\right)\right)$ as

$\begin{array}{l}\left(\mathcal{S}u\right)\left(t\right)={\displaystyle {\int }_0^1}{K}_s\left(t,s\right)u\left(s\right)\text{d}s\text{with}{K}_s\left(t,s\right)=\frac{\left|x^{\prime }\left(s\right)\right|}{2\pi }\ln \frac{1}{\left|x\left(t\right)-x\left(s\right)\right|},\\ \left(\mathcal{D}u\right)\left(t\right)={\displaystyle {\int }_0^1}{K}_c\left(t,s\right)u\left(s\right)\text{d}s\text{with}{K}_c\left(t,s\right)=\{\begin{array}{ll}\frac{1}{2\pi }\frac{x^{\prime }{\left(s\right)}^{\perp }\cdot \left(x\left(t\right)-x\left(s\right)\right)}{{\left|x\left(t\right)-x\left(s\right)\right|}^2},\hfill & t\ne s,\hfill \\ \frac{1}{4\pi }\frac{x^{\prime }{\left(t\right)}^{\perp }\cdot x^{″}\left(t\right)}{{\left|x^{\prime }\left(t\right)\right|}^2},\hfill & t=s,\hfill \end{array}\\ \left({\mathcal{D}}^{\prime }u\right)\left(t\right)={\displaystyle {\int }_0^1}{K^{\prime }}_c\left(t,s\right)u\left(s\right)\text{d}s\text{with}{K^{\prime }}_c\left(t,s\right)=K_c\left(s,t\right)\frac{\left|x^{\prime }\left(s\right)\right|}{\left|x^{\prime }\left(t\right)\right|}\end{array}$

for $0\le t\mathrm{,}s\le 1$, where we denote ${\left(x_1,x_2\right)}^{\perp }=\left(x_2,-x_1\right)$. The kernel $K_s$ is weakly singular while $K_c$ and ${K^{\prime }}_c$ are continuous. The singularity in $K_s$ can be rearranged as

$\ln \left|x\left(t\right)-x\left(s\right)\right|=\ln \left(2\left|\sin \left(\pi \left(t-s\right)\right)\right|\right)+K_0\left(t,s\right)$

with continuous kernel

$K_0\left(t,s\right)=\{\begin{array}{ll}\ln \frac{\left|x\left(t\right)-x\left(s\right)\right|}{2\left|\sin \left(\pi \left(t-s\right)\right)\right|},\hfill & t\ne s\hfill \\ \ln \frac{\left|x^{\prime }\left(t\right)\right|}{2\pi },\hfill & t=s\hfill \end{array}$

so that integrals involving this singularity can be dealt with by exact integration, as we employ Nyström’s method with trigonometric interpolation on regular grids (see e.g. [12], Chapter 12).

In our examples, we take $\Omega$ as the elliptic region for the sake of simplicity:

$\Omega =\left\{\left(x_1,x_2\right):x_1^2/a^2+x_2^2/b^2<1\right\}\text{with}\left(a,b\right)=\left(1,0.4\right).$

We should note that this choice of the aspect ratio of the ellipse gives rise to a more severely ill-posed case of the Robin inverse problem than most cases in other studies, yet as observed consistently in our numerical experiments and demonstrated in examples below, the inclusion of nonnegativity constraints in the iterative schemes does help to combat the ill-posedness effectively and produce markedly better results for p.

Standard parametrization of the ellipse is used:

$x=x\left(t\right)=\left(a\cos \left(2\pi t\right),b\sin \left(2\pi t\right)\right),0\le t\le 1.$

The two segments ${\Gamma }_1$ (where p is defined) and ${\Gamma }_0$ (where $u_0$ is given) are chosen as

${\Gamma }_1=\left\{x\left(t\right):t\in \left[0.1,0.4\right]\right\}\text{and}{\Gamma }_0=\left\{x\left(t\right):t\in \left[0.6,0.9\right]\right\}.$

The input function g is set as

$g\left(t\right)=\{\begin{array}{ll}1,\hfill & \text{when}t\in \left[0.5,0.6\right]\hfill \\ 0,\hfill & \text{elsewhere}\text{in}\left[0,1\right],\hfill \end{array}$

and the true profile for p as

$p\left(t\right)=\{\begin{array}{ll}{\sin }^4\left(\frac{t-0.1}{0.3}\pi \right),\hfill & \text{when}t\in \left[0.1,0.4\right]\hfill \\ 0,\hfill & \text{elsewhere}\text{in}\left[0,1\right].\hfill \end{array}$

Discretization grid size is set to $h=1/400$ for discretization of the integral operators and the systems $\mathcal{A}w=f$ in (3.4) and $\mathcal{B}\left(z\right)=f$ in (3.9). To generate the synthetic data $u_0$ on ${\Gamma }_0$, we solve the discretized alternative formulation (2.3) - (2.4), and add to it uniformly distributed random noise of noise level $\delta =0.01$, relative to the L² norm of $u_0$.

In both iterative methods, we stop the iteration when the L² norm of the Gauss-Newton step $d^{\left(k\right)}$ from (3.7) or (3.11) is less than 10⁻¹⁰.

Our selection of the regularization parameter $\alpha$ is guided by the idea of the normalized cumulative periodogram (NCP) method (see e.g. [15]), which is based on monitoring the noise pattern in the data-fitting residual ${\mathcal{R}}_{0}u-u_0$ in (3.4) or (3.9). If the residual is dominated by white noise, then $\alpha$ is a suitable choice; if it is dominated by high frequency noise, the current parameter $\alpha$ shall be increased, and if it is dominated by a low frequency signal, then $\alpha$ shall be reduced to a smaller value ( [16]). See also ( [10], Chapter 7) for other methods of selecting a proper regularization parameter value for illposed problems.

The performance of the algorithms is less sensitive to the choice of $\beta$, and they work very similarly for values of $\beta$ from the range between 10⁻¹ to 10². We set $\beta =1$ in our experiments.

We first experiment using the linear least-squares method (3.7) - (3.8), and the result of one such example is presented in Figure 1. We select $\alpha =5.0\times {10}^{-9}$. Starting with $w^{\left(0\right)}=0$, it typically takes only a few iterations (less than 10) to satisfy the stopping criterion throughout our experiments. For comparison, we also present the result without the penalty term (when $\beta =0$ ), which is direct and does not require iterations, but followed by converting from $\left(u\mathrm{,}v\right)$ to p

using $p=u^{+}{v}^{+}/\left(\epsilon +{\left(u^{+}\right)}^2\right)$, as in [11], where $f^{+}=\max \left(f,0\right)$ and $\epsilon >0$ as

the machine epsilon. As seen, the few iteration needed for the nonnegativity constraints indeed have helped to improve the overall quality of the recovered coefficient p.

We then test on the nonlinear least-squares method (3.11) - (3.12), with the result of one such example presented in Figure 2. We choose $\alpha =5.0\times {10}^{-10}$ in

Figure 1. Recovered p by the linear least-squares method with nonnegativity constraints (solid), and by the direct method with truncation (dotted), in comparison with the true profile (dashed).

Figure 2. Recovered p by the nonlinear least-squares method, with nonnegativity constraints (solid) and without (dotted), in comparison with the true profile (dashed).

this case. Typically it takes about 10 iterations from initial $z^{\left(0\right)}=0$ to satisfy the stopping criterion in our experiments. For comparison, we also carry out parallel calculations without the penalty term (by setting $\beta =0$ ), and as expected, the resulting solution assumes some negative values due to the absence of nonnegativity constraints. As can be seen, adding the penalty term not only ensures the nonnegativity of the recovered coefficient p, but also improves its quality.

Throughout our numerical experiments, we have consistently observed the effectiveness of this approach in dealing with the ill-posedness of the Robin inverse problem and in improving recovery results, which is achieved within the same iterative framework and without additional computational cost.

We conclude with some final remarks: 1) It is possible to fine-tune these least-squares methods by introducing different weights between the model equation residual and the data fitting residual in (3.4) or (3.9), but at the expense of one more parameter (the weight ratio) to adjust, as studied in [14]. Even in this setting, the penalty strategy proposed here remains readily applicable, and the penalty parameter $\beta$ is generally less sensitive and easier to determine than the weight ratio and regularization parameter; 2) In the more common approach where the dependence of u on p is determined exactly by solving the model equation (1.1) or (2.2) and only the data fitting residual is to be minimized, it is still possible to include the penalty strategy proposed here to better facilitate the solution process for maintaining nonnegativity of p, and it is reasonable to expect the improvement of quality of the recovered p in more ill-posed situations, as we have experienced and reported here.

Acknowledgements

This work was supported by the key research projects of general universities in Guangdong Province (2019KZDXM034), and the basic research and applied basic research projects in Guangdong Province (Projects of Guangdong-Hong Kong-Macao Center for Applied Mathematics) (2020B1515310018).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Loh, W.H. (1987) Modeling and Measurement of Contact Resistances. Technical Report No. G830-1, Stanford Electronics Labs, Stanford.
[2]	Fang, W. and Cumberbatch, E. (1992) Inverse Problems for MOSFET Contact Resistivity. SIAM Journal on Applied Mathematics, 52, 699-709. https://doi.org/10.1137/0152039
[3]	Fang, W. and Lu, M. (2004) A Fast Wavelet Collocation Method for an Inverse Boundary Value Problem. International Journal for Numerical Methods in Engineering, 59, 1563-1585. https://doi.org/10.1002/nme.928
[4]	Kaup, P. G. and Santosa, F. (1995) Nondestructive Evaluation of Corrosion Damage Using Electrostatic Measurements. Journal of Nondestructive Evaluation, 14, 127-136. https://doi.org/10.1007/BF01183118
[5]	Inglese, G. (1997) An Inverse Problem in Corrosion Detection. Inverse Problems, 13, 977-994. https://doi.org/10.1088/0266-5611/13/4/006
[6]	Chaabane, S. and Jaoua, M. (1999) Identification of Robin Coefficients by Means of Boundary Measurements. Inverse Problems, 15, 1425-1438. https://doi.org/10.1088/0266-5611/15/6/303
[7]	Jin, B. (2007) Conjugate Gradient Method for the Robin Inverse Problem Associated with the Laplace Equation. International Journal for Numerical Methods in Engineering, 71, 433-453. https://doi.org/10.1002/nme.1949
[8]	Cakoni, F. and Kress, R. (2007) Integral Equations for Inverse Problems in Corrosion Detection from Partial Cauchy Data. Inverse Problems and Imaging, 1, 229-245. https://doi.org/10.3934/ipi.2007.1.229
[9]	Lin, F.-R. and Fang, W. (2005) A Linear Integral Equation Approach to the Robin Inverse Problem. Inverse Problems, 21, 1757-1772. https://doi.org/10.1088/0266-5611/21/5/015
[10]	Vogel, C. (2002) Frontiers in Applied Mathematics: Computational Methods for Inverse Problems. Society for Industrial and Applied Mathematics, Philadelphia. https://doi.org/10.1137/1.9780898717570
[11]	Fang, W. and Zeng, S. (2009) A Direct Solution for the Robin Inverse Problem. Journal of Integral Equations and Applications, 21, 545-557. https://doi.org/10.1216/JIE-2009-21-4-545
[12]	Kress, R. (1999) Linear Integral Equations. 2nd Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-0559-3
[13]	McLean, W. (2000) Strongly Elliptic Systems and Integral Equations. Cambridge University Press, Cambridge.
[14]	Fang, W., Lin, F.-R. and Ma, Y. (2019) Fast algorithms for Boundary Integral Equations on Elliptic Domains and Related Inverse Problems. East Asian Journal on Applied Mathematics, 9, 485-505. https://doi.org/10.4208/eajam.130917.170818
[15]	Hansen, P.C., Kilmer, M.E. and Kjeldsen, R.H. (2006) Exploiting Residual Information in the Parameter Choice for Discrete Ill-Posed Problems. BIT Numerical Mathematics, 46, 41-59. https://doi.org/10.1007/s10543-006-0042-7
[16]	Fang, X., Lin, F.-R. and Wang, C. (2017) Estimation of a regularisation Parameter for a Robin Inverse Problem. East Asian Journal on Applied Mathematics, 7, 325-342. https://doi.org/10.4208/eajam.150216.260117a