On the Non-Linear Integral Equation Approach for an Inverse Boundary Value Problem for the Heat Equation

TL;DR

This paper presents a novel non-linear integral equation approach for solving an inverse boundary problem in heat conduction, utilizing Laguerre transforms, boundary integral equations, and regularization techniques for accurate reconstructions.

Contribution

It introduces a new method combining Laguerre transforms and non-linear boundary integral equations for inverse heat problems, with a developed iterative algorithm and stability analysis.

Findings

The method achieves accurate boundary reconstructions.

The approach demonstrates stability through Tikhonov regularization.

Numerical results confirm the effectiveness of the proposed technique.

Abstract

We consider the inverse problem of reconstructing the interior boundary curve of a doubly connected domain from the knowledge of the temperature and the thermal flux on the exterior boundary curve. The use of the Laguerre transform in time leads to a sequence of stationary inverse problems. Then, the application of the modified single-layer ansatz, reduces the problem to a sequence of systems of non-linear boundary integral equations. An iterative algorithm is developed for the numerical solution of the obtained integral equations. We find the Fr\'echet derivative of the corresponding integral operator and we show the unique solvability of the linearized equation. Full discretization is realized by a trigonometric quadrature method. Due to the inherited ill-possedness of the derived system of linear equations we apply the Tikhonov regularization. The numerical results show that the proposed method produces accurate and stable reconstructions.

Equations

$$\frac{1}{\alpha }\frac{\partial u}{\partial t}=\Delta u,\text{mbox}in\Omega \times (0,\infty )$$

$$u(\cdot ,0)=0,\text{mbox}in\Omega$$

$$u=0,\text{mbox}on{\Gamma }_1\times [0,\infty ),\frac{\partial u}{\partial \nu }=g,\text{mbox}on{\Gamma }_2\times [0,\infty ).$$

$$u=f,\text{mbox}on{\Gamma }_2\times (0,\infty ).$$

$$u(x,t)=\kappa \sum _{n=0}^{\infty }{u}_n(x)L_n(\kappa t),$$

$$u_n(x):={\int }_0^{\infty }{e}^{-\kappa t}{L}_n(\kappa t)u(x,t)dt,n=0,1,2,\dots$$

$$L_n(t)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)\frac{(-t{)}^k}{k!},$$

$$(n+1)L_{n+1}(t)$$

$$L_{n+1}^{\prime }(t)$$

$$\Delta u_n-{\gamma }^2{u}_n=\beta \sum _{m=0}^{n-1}{u}_m,\text{mbox}in\Omega$$

$$u_n=0,\text{mbox}on{\Gamma }_1\text{and}\frac{\partial u_n}{\partial \nu }=g_n,\text{mbox}on{\Gamma }_2.$$

$$g_n(x):={\int }_0^{\infty }{e}^{-\kappa t}{L}_n(\kappa t)g(x,t)dt,n=0,1,2,\dots$$

$$u_n=f_n,\text{mbox}on{\Gamma }_2$$

$$f_n(x):={\int }_0^{\infty }{e}^{-\kappa t}{L}_n(\kappa t)f(x,t)dt,n=0,1,2,\dots .$$

$$I_0(z)=\sum _{n=0}^{\infty }\frac{1}{(n!{)}^2}{\left(\frac{z}{2}\right)}^{2n},I_1(z)=\sum _{n=0}^{\infty }\frac{1}{n!(n+1)!}{\left(\frac{z}{2}\right)}^{2n+1},$$

$$K_0(z)$$

$$K_1(z)$$

$$\psi (n)=\sum _{m=1}^n\frac{1}{m},n=1,2,\dots$$

$$v_n(r)=\sum _{k=0}^{\left[\frac{n}{2}\right]}{a}_{n,2k}{r}^{2k},w_n(r)=\sum _{k=0}^{\left[\frac{n-1}{2}\right]}{a}_{n,2k+1}{r}^{2k+1}$$

$$a_{n,0}$$

$$a_{n,n}$$

$$a_{n,k}$$

$${\Phi }_n(x,y):=K_0(\gamma |x-y|)v_n(|x-y|)+K_1(\gamma |x-y|)w_n(|x-y|),x\ne y,$$

$$u_n(x)=\frac{1}{\pi }\sum _{\ell =1}^2\sum _{m=0}^n{\int }_{{\Gamma }_{\ell }}{\varphi }_m^{\ell }(y){\Phi }_{n-m}(x,y)ds(y),x\in \Omega$$

$$\frac{1}{\pi }\sum _{\ell =1}^2{\int }_{{\Gamma }_{\ell }}{\varphi }_n^{\ell }(y){\Phi }_0(x,y)ds(y)$$

$${\varphi }_n^2(x)+\frac{1}{\pi }\sum _{\ell =1}^2{\int }_{{\Gamma }_{\ell }}{\varphi }_n^{\ell }(y)\frac{\partial {\Phi }_0}{\partial n(x)}(x,y)ds(y)$$

$$\frac{1}{\pi }\sum _{\ell =1}^2{\int }_{{\Gamma }_{\ell }}{\varphi }_n^{\ell }(y){\Phi }_0(x,y)ds(y)$$

$$F_{1,n}(x)$$

$$G_n(x)$$

$$F_{2,n}(x)$$

$${\Gamma }_{\ell }=\{x_{\ell }(s)=(x_{1\ell }(s),x_{2\ell }(s)),s\in [0,2\pi ]\},\ell =1,2$$

$${\phi }_n^{\ell }(s):={\varphi }_n^{\ell }(x_{\ell }(s))|x_{\ell }^{\prime }(s)|.$$

$$\frac{1}{2\pi }\sum _{\ell =1}^2{\int }_0^{2\pi }{\phi }_n^{\ell }(\sigma )H_0^{1,\ell }(s,\sigma )d\sigma$$

$$\frac{{\phi }_n^2(s)}{|x_2^{\prime }(s)|}+\frac{1}{2\pi }\sum _{\ell =1}^2{\int }_0^{2\pi }{\phi }_n^{\ell }(\sigma )Q_0^{2,\ell }(s,\sigma )d\sigma$$

$$\frac{1}{2\pi }\sum _{\ell =1}^2{\int }_0^{2\pi }{\phi }_n^{\ell }(\sigma )H_0^{2,\ell }(s,\sigma )d\sigma$$

Full text

Keywords: doubly connected domain; boundary reconstruction; Laguerre transform; modified single layer potentials; non-linear boundary integral equations; trigonometric quadrature method; Newton method; Tikhonov regularization.

1 Introduction

Integral equation methods are among the most frequently used and effective methods for solving problems in applied sciences. The main advantage is that they are methods for dimension reduction and instead of solving a problem in an unbounded region, we obtain the solution from a boundary integral equation. This has many applications to a wide class of direct boundary value problems but also to various inverse problems.

In this work, we deal with the problem of reconstructing the inner boundary of a medium. There are different ways to solve this problem. Once the inverse problem is written as a non-linear operator equation, then its approximate solution can be found by Newton’s method (see i.e. [8, 10, 11]). The inverse problem can also be directly reduced to a system of non-linear integral equations. This is the reciprocity gap approach, based on the Green’s integral theorem. In [14, 15, 18], it was used for the case of the Laplace and Helmholtz equations. Another way to reduce the problem of reconstructing the boundary curve of a medium to a system of non-linear integral equations is to use potentials for the integral representation of the solution [1, 3, 6, 12, 16]. Representing the unknown function in a form of, for example, a single-layer potential, the inverse problem is reduced to a system of non-linear equations with respect to the unknown densities and the unknown boundary of the domain. Then, we may apply an iterative scheme.

In the case of non-linear inverse problem for the non-stationary heat conductivity equation, an approach with the use of the indirect method of integral equations for the solution of the corresponding non-linear operator equation was developed in [9, 10]. In this paper, we want to use our experience with solving non-stationary direct problems by a combination of the Laguerre transform with integral equations [4, 7, 8]. The application of the Laguerre transform with respect to the time variable, results to the semi-discretization of the given non-linear inverse parabolic problem. This allows us to obtain a sequence of stationary inverse problems, which is then further reduced by the method of integral potentials to a sequence of systems of non-linear integral equations.

The paper is organized as follows. In the following, we present the governing equations and we formulate the corresponding inverse problem. In section 2, we apply the Laguerre transform to the inverse boundary problem with respect to the time variable and we reduce it to a sequence of stationary boundary problems. Then, using the modified single-layer ansatz, we rewrite it equivalently as a sequence of non-linear boundary integral equations. Here, the proposed iterative scheme for its solution is also described. The section 3 contains the numerical solution of system of the integral equations. Firstly, we consider the well-posed system of the “field” equations and we apply the Nyström method for its discretization. Then, we linearize the ill-posed “data” equation and apply the collocation method for its discretization. The Tikhonov regularization is used for solving the derived ill-conditioned linear system. Finally, in section 4 we provide numerical examples where the reconstruction of the unknown boundary is feasible for different setups.

Let $\Omega$ be a doubly connected domain in ${\rm I\!R}^{2}$ with smooth boundary $\Gamma$ of class $C^{2}.$ We assume that $\Gamma$ consists of two disjoints curves $\Gamma_{1}$ and $\Gamma_{2}$, meaning $\Gamma=\Gamma_{1}\cup\Gamma_{2},$ with $\Gamma_{1}\cap\Gamma_{2}=\emptyset,$ such that $\Gamma_{1}$ is contained in the interior of $\Gamma_{2}.$

We consider the following initial boundary value problem for the heat equation

[TABLE]

subject to the homogeneous initial condition

[TABLE]

and the boundary conditions

[TABLE]

Here $\alpha$ is the thermal diffusivity, $\nu$ denotes the outward unit normal to $\Gamma$ and $g$ is a given and sufficiently smooth function. Existence and uniqueness of classical (see [13, 17]) and of weak solutions (see [19, 20]) of this mixed initial boundary value problem is well established.

In this work we consider the inverse problem of determining the interior boundary curve $\Gamma_{1}$ from the knowledge of the thermal flux $g$ and of the temperature

[TABLE]

2 Laguerre transform and Boundary Integral Equations

We solve the system (1.1) – (1.4) in two steps. First, we represent the solution $u$ as a (scaled) Fourier expansion with respect to the Laguerre polynomials resulting to a sequence of stationary mixed problems. Then, we reduce the problem to a system of boundary integral equations.

We consider the expansion

[TABLE]

where

[TABLE]

for $\kappa>0,$ and the Laguerre polynomial $L_{n}$ of order $n,$ defined by

[TABLE]

and satisfying the recurrence relations

[TABLE]

for $n=0,1,2,\ldots.$ Using these formulas, we can show (see [10]) that the Fourier–Laguerre coefficients $u_{n}$, satisfy the following sequence of mixed problems

[TABLE]

with boundary conditions

[TABLE]

Here,

[TABLE]

are the Laguerre–Fourier coefficients of the given function $g$ and $\beta=\kappa/\alpha,\,\gamma^{2}=\beta.$ The additional condition (1.4) results to

[TABLE]

with

[TABLE]

In order to apply the non-linear integral equation method we need the fundamental solution of the equations (2.1). We consider the modified Bessel functions

[TABLE]

and the modified Hankel functions

[TABLE]

of order zero and one, respectively. Here, we set $\psi(0)=0$ and

[TABLE]

and $C=0.57721\ldots$ denotes the Euler constant (see [2]).

We define the polynomials $v_{n}$ and $w_{n}$ by

[TABLE]

with the convention $w_{0}(r)=0.$ The coefficients are given by the relations

[TABLE]

for $n=1,2,\ldots$.

Then, the sequence of functions

[TABLE]

satisfies (2.1) with respect to $x$ in ${\rm I\!R}^{2}\setminus\{y\}$ for $n=0,1,2,\ldots,$ that is, the function $\Phi_{n}$ provides a fundamental solution (see [4]). Now we can represent the solutions $u_{n}$ of the problem (2.1) – (2.3) using the following single layer potential form

[TABLE]

with the unknown densities $\phi_{m}^{1}$ and $\phi_{m}^{2}$, $m=0,1,\ldots$, defined on the boundary curves $\Gamma_{1}$ and $\Gamma_{2}$, respectively, and $\Phi_{n}$ is given by (2.6).

We let $x$ tend to the boundary and using the boundary conditions and the standard jump relations we get the following system

[TABLE]

with right-hand sides

[TABLE]

This is a system of three equations for the three unknowns: the two densities $\phi^{1}_{n},\,\phi^{2}_{n}$ and the boundary $\Gamma_{1}.$ The operators are singular, linear on the densities but act non-linearly on the boundary curve. We will consider the Frechét derivative of the integral operators for linearizing it.

There could be many ways to solve this system of integral equations. Here, motivated by [16, 18], we propose the following iterative scheme which avoids the full linearization of the system and the computation of the Fréchet derivative of the normal derivative operator, which is more involved compared to the single layer operator. We refer to [1, 3, 12] for some recent applications of this scheme to different regimes.

Before presenting the iterative method, we consider the parametrization of the system (2.8) – (2.10). We assume the following parametric representation of the boundary

[TABLE]

and we define

[TABLE]

Then, the system (2.8) – (2.10) takes the form

[TABLE]

for $n=0,\ldots,N$, $N\in{\rm I\!N}$ and right-hand sides

[TABLE]

The kernels are given by

[TABLE]

for $s\neq\sigma,\,\,k,\ell=1,2$, and $n=0,\ldots,N$. The functions $\Phi_{n}$ are defined in (2.6).

We solve the above system of equations using the iterative scheme:

Step 1

Given an initial approximation of $\Gamma_{1}$, we solve the sequence of well-posed systems of integral equations (2.11) and (2.12) for $\phi^{1}_{n},\,\phi^{2}_{n},\,\,n=0,...,N.$

Step 2

Keeping now the densities fixed, we linearize the ill-posed integral equation (2.13) resulting to

[TABLE]

We solve the $n=0,...,N$ equations for the parametrization $\chi$ of the perturbed $\Gamma_{1},$ and we update as $x_{1}+\chi.$

Equation (2.15) contains the Fréchet derivative $\mathcal{D}_{n}$ of the integral operator with kernel $H^{2,1}_{n}$ with respect to $x_{1}$ as a linear operator on $\chi.$ The analytic form of this operator is given in the next section. In (2.15) we have $N$ equations to be solved for one function $\chi.$ We present later two different ways for solving this equation.

3 Numerical Implementation

In this section we consider the numerical implementation of the iterative scheme presented above. We refer to (2.11) – (2.12) as the “field” system and to (2.13) as the “data” equation.

3.1 Numerical solution of the “field” system

The singular kernels $H^{\ell,\ell}_{n}$ appearing in (2.11) and in (2.13), are analyzed using the explicit expression (2.6) for the elements of the fundamental sequence together with the expansions (2.4) – (2.5). Then, we get the expression

[TABLE]

After lengthy but straightforward calculations, we treat the logarithmic singularity, using the following representation

[TABLE]

where

[TABLE]

and

[TABLE]

The diagonal terms are given by

[TABLE]

For the representation of the kernels $Q^{k,\ell}_{n}$, we introduce the function

[TABLE]

and the polynomials

[TABLE]

Then, the kernels $Q^{k,\ell}_{n}$ admit the form

[TABLE]

for $s\neq\sigma$ and $n=0,1,\ldots,N$.

The kernels $Q^{\ell,\ell}_{n}$ have logarithmic singularity. As in the case of the kernels $H^{\ell,\ell}_{n}$ performing similar calculations, we derive the decomposition

[TABLE]

where

[TABLE]

and

[TABLE]

The diagonal elements are now given by

[TABLE]

Now we have the explicit representations of the singular kernels, and we can apply the following standard quadrature rules [17] for the numerical discretization

[TABLE]

for the mesh points

[TABLE]

using the weight function

[TABLE]

Thus, we approximate the solution of the “field” system (2.11) – (2.12) by collocating the integral equations at the nodal points $\{s_{k}\}$ leading to the sequence of linear systems

[TABLE]

for $i=0,\ldots,2M-1$, with the right-hand sides

[TABLE]

and

[TABLE]

where we used the abbreviation $\varphi_{n,j}^{\ell}\approx\varphi_{n}^{\ell}(s_{j})$, $\ell=1,2$, $n=0,\ldots,N$, $j=0,\ldots,2M-1$.

3.2 Numerical solution of the “data” equation

As described in our algorithm, we search for the correction of $\Gamma_{1}$ by solving the “data” equation (2.13), assuming that we know the densities $\varphi^{\ell}_{n}$, $\ell=1,2$, $n=0,...,N$.

For simplicity, we consider starlike interior curve, meaning we assume parametrization in polar coordinates of the form

[TABLE]

where $r:{\rm I\!R}\to(0,\infty)$ is a $2\pi$ periodic function representing the radial distance from the origin. However, we wish to stress that the following analysis is also applicable to other boundaries.

The linearized equation (2.15) admits the following parametric form

[TABLE]

where $q$ is the radial function of the perturbed boundary. The Fréchet derivative $\mathcal{D}_{n}$ has the explicit form

[TABLE]

with kernel

[TABLE]

where

[TABLE]

Theorem 3.1

The Fréchet derivative operator $\mathcal{D}_{n}[\varphi,r;q]$ is injective at the exact solution.

*Proof. *We just have to show that if $q$ solves

[TABLE]

then $q=0.$ We follow the ideas of [5] and we set

[TABLE]

where $\zeta(x_{1}(s))=q(s)(\cos s,\sin s)$. Clearly $V_{n}$ satisfies the sequence (2.1). We define $V_{n}^{\pm}(x):=\lim_{h\rightarrow 0^{+}}V_{n}(x\pm h\nu).$ Then, by assumption $V_{n}^{+}|_{\Gamma_{1}}=0$. The perturbed interior curve can be represented as follows [5, 15]

[TABLE]

for small perturbations and a given function $\tilde{q}.$

Then we can rewrite the functions $V_{n}$ in the form

[TABLE]

By the properties of the double layer potentials $V_{n}$ [8], we extend it continuously to $\Gamma_{1}$ as

[TABLE]

By the uniqueness of the exterior and interior Dirichlet problems we have

[TABLE]

The functions $u_{n}$ given by (2.7) solve the Dirichlet problem in the interior of $\Gamma_{1}$ with homogeneous boundary conditions. Then by the unique solvability, the functions $u_{n}$ have to vanish in the interior of $\Gamma_{1}$ and hence $\tfrac{\partial u_{n}^{-}}{\partial\nu}=0$ on $\Gamma_{1}$. The jump relations imply that $\frac{\partial u_{n}^{+}}{\partial\nu}|_{\Gamma_{1}}=\sum_{m=0}^{n}\varphi_{m}^{1}(s)$. Employing Holmgren’s uniqueness theorem to the sequence of Helmholtz equations, one can show that the Cauchy data $(u_{n}^{+},\,\frac{\partial u_{n}^{+}}{\partial\nu})$ cannot be identically zero on an open subset of $\Gamma_{1}$. Thus, in view of (3.3), we conclude that $\tilde{q}=0$ and therefore $q=0.$   $\Box$

We apply the quadrature rules, analyzed in the previous section, in (3.2) and then a collocation method approximating the function $q$ by a trigonometric polynomial of the form

[TABLE]

with

[TABLE]

We substitute (3.4) in (3.2) and at the nodal points $\{s_{i}\}$ we obtain the following linear system

[TABLE]

where \mbox{\boldmath{q}}=(q_{0},...,q_{2J})^{\top}\in{\rm I\!R}^{2J+1}, and \mbox{\boldmath{A}}\in{\rm I\!R}^{(2M)\times(2J+1)},\,\mbox{\boldmath{b}}\in{\rm I\!R}^{2M} given by

[TABLE]

We solve the linear system (3.5) either for $n=N,$ or using all available information, meaning

[TABLE]

We refer to these cases as the “final-step” and the “multi-step” system, respectively.

Due to the ill-posedness of (3.5) and its overdetermination we apply the least-squares method with Tikhonov regularization. Then, we solve

[TABLE]

with the regularization parameter $\lambda>0,$ to be chosen by trail and error.

4 Numerical Results

We present numerical examples for different boundary curves. We consider the following three cases:

Example 1

The exterior boundary curve $\Gamma_{2}$ is a circle with center $(0,0)$ and radius $1,$ and the interior boundary curve $\Gamma_{1}$ (to be reconstructed) is peanut-shaped with radial function

[TABLE]

Example 2

The exterior boundary curve $\Gamma_{2}$ is a rounded rectangle with radial function

[TABLE]

and $\Gamma_{1}$ is a apple-shaped boundary with radial function

[TABLE]

Example 3

Both boundary curves are kite-shaped with parametrizations

[TABLE]

and

[TABLE]

We generate the simulated Cauchy data by solving the sequence

[TABLE]

with boundary conditions

[TABLE]

for given boundary functions $f_{\ell,n},\,\ell=1,2.$ To avoid an inverse crime, we consider double amount of nodal points compared to the inverse problem and we compute the thermal flux $g$ by tending the solution to the exterior boundary and considering the jump relations. Then, we add error to the Cauchy data with respect to the $L^{2}$ norm

[TABLE]

for given noise level $\delta,$ and normally distributed random variables $u,v\in{\rm I\!R}.$ At every step, we update the regularization parameter using

[TABLE]

The boundary functions are chosen as

[TABLE]

In all examples, the initial guess is a circle with center $(0,0)$ and radius $r_{0},$ we set $\alpha=1$ and $M=64.$ We observed that the additional information provided by the “multi-step” system produces more stable reconstructions for noisy data and requires less iterations for exact data. Thus, in the following examples we consider the “final-step” method for exact data and the “multi-step” method for noisy data. In the first example, we set $\kappa=1$ and $N=10.$ The results are presented in Figure 1, for $J=5$ coefficients, using $r_{0}=0.4$ as initial radius. The reconstruction for noise-free data are obtained using as initial regularization parameter $\lambda_{0}=0.0001$ and 24 iterations. For data with $3\%$ noise, we used $\lambda_{0}=0.001$ and 14 iterations.

In Figure 2, we see the reconstructions for the parametrized boundary curves of Example 2. We keep all the parameters the same as in the first example. We mark here that the results are again satisfactory independently of the boundary parametrizations. The regularization parameters are also kept fixed and the reconstructions are after 20 iterations (exact data) and after 9 iterations (data with $3\%$ noise).

In the third example, we test the performance of the algorithm for more irregular boundary curves. Here, we set $J=7$ and the initial radius is $r_{0}=0.5.$ The results presented in Figure 3 are after 13 iterations for exact data and after 9 iterations for noisy data.

5 Conclusion

We developed a non-linear integral equations approach for the inverse parabolic problem related to the reconstruction of a part of the boundary curve from the knowledge of the Cauchy data on the other part of the boundary. Our strategy consisted on the consecutive dimension reduction. Firstly we considered the semi-discretization by Laguerre transform in time for the given three-dimensional problem. It resulted to a sequence of stationary inverse boundary problems for the Helmholtz equation. With the help of the modified single layer potentials representation, these problems were reduced to a sequence of one-dimensional non-linear boundary integral equations. Then, a Newton-type iteration method was applied. We solved the well-posed system of linear integral equations by the Nyström method and the ill-posed linear integral equation by the collocation method with Tikhonov regularization, at every iteration step. In general, our approach can be applied without significant changes to the inverse boundary problems for the hyperbolic equation and can also be extended to the case of three-dimensional domains.