Identifying an unknown source in the Poisson equation by a wavelet dual least square method

Abstract

This paper deals with an inverse problem of identifying an unknown source which depends only on one variable in two-dimensional Poisson equation, with the aid of an extra measurement at an internal point. This problem is ill-posed, we proposed a regularization strategy, a wavelet dual least square method, to analyze the stability of the problem. Meanwhile, a numerical experiment is devised to verify the validity of the method.

MSC:35R40, 65J20.

1 Introduction

Consider the following inverse problem: find a pair of functions $(u(x,y),f(x))$ satisfying

$$\{\begin{array}{cc}-u_{xx}-u_{yy}=f(x),\hfill & -\mathrm{\infty }<x<\mathrm{\infty },y>0,\hfill \\ u(x,0)=0,\hfill & -\mathrm{\infty }<x<\mathrm{\infty },\hfill \\ u(x,y){|}_{y\to \mathrm{\infty }}\text{ bounded},\hfill & -\mathrm{\infty }<x<\mathrm{\infty },\hfill \\ u(x,1)=g(x),\hfill & -\mathrm{\infty }<x<\mathrm{\infty },\hfill \end{array}$$

(1.1)

where $f(x)$ is the unknown source depending only on one spatial variable and $u(x,1)=g(x)$ is the supplementary condition. In practical applications, the input data $g(x)$ can only be measured. There will be measured data function $g_{\delta }(x)$ which is merely in $L^2(R)$ and satisfies

$${\parallel g-g_{\delta }\parallel }_{L^2(R)}\le \delta ,$$

(1.2)

where the constant $\delta >0$ represents a noise level of input data. This problem is called the inverse problem of unknown source identification.

Inverse source identification problems are important in many branches of engineering sciences such as crack determination [1, 2], heat source determination [3], heat conduction problems [4, 5], electromagnetic theory [6]. This kind of problem arises in many important applications in practice, e.g., with the development of society and economics, groundwater pollution has become a serious threat to the environment. The government has to take some measures to prevent the groundwater from further contaminations. But the cost of cleanup for polluted aquifers is staggering, and in many cases it is hard to identify which companies are responsible for the contamination due to the lack of tools to discover the pollution sources. So, it is necessary to try to give more concrete information of the characteristics (location, magnitude, and duration of activity) of specific groundwater pollution sources (see [7]). As we know, most attempts at quantifying contaminant transport rely on mathematical methods. Since the data cannot be measured by direct ways in many cases, we are always encountering inverse problems of deciding unknown sources and aquifer parameters.

The investigation of the traditional inverse potential problem can be found in [8, 9]. The studies of such problems give a complete analysis of experimental data. In general, a full source f in (1.1) is not solely attainable from boundary measurements. The inverse source identification problem becomes solvable if some a priori knowledge is assumed. For instance, when one of the products in the separation of variables is known [7, 10], or the base area of a cylindrical source is known [10], or a non-separable type is in the form of a moving front [7], the boundary data g can then uniquely determine the unknown sources f. Furthermore, when both u and f are relatively smooth, some standard regularization techniques can be employed (see [11] for a more detailed overview).

Wavelet regularization methods have been studied for solving various types of inverse problems in the heat equation [12–16]. Eldén [12] and Regińska [13, 14], Xiong [15] used the wavelet Galerkin method and the wavelet method to approximate the sideways heat equation by Meyer wavelets, and Xiong [16] used the wavelet dual least squares method to approximate the BHCP by Shannon wavelets. In this work, by using Meyer wavelets, we obtain an explicit error estimate of Hölder type between the unknown source term and its approximation. Moreover, according to the general theory of regularization, we conclude that our estimate is order optimal.

In general, for an ill-posed problem, the convergence rates of the regularization solution can only be given under a priori assumptions on the exact data. We will formulate such an a priori assumption in terms of an exact solution $f(x)$ by considering

$${\parallel f\parallel }_p\le E,p>0,$$

(1.3)

where the ${\parallel \cdot \parallel }_p$ denotes the Sobolev space $H^p(\mathbb{R})$-norm defined by

$${\parallel f\parallel }_p:={\left({\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{(1+{\xi }^2)}^p{|\hat{f}(\xi )|}^{2}d\xi \right)}^{1/2}.$$

(1.4)

In order to formulate problem (1.1) in terms of an operator equation in the space $X=L^2(R)$, let A be the operator on X defined as follows:

$$Af(x)=g(x).$$

(1.5)

Let

$$\hat{f}(\xi )=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)e^{-i\xi x}dx,\xi \in R$$

(1.6)

be the Fourier transform of the function $f(x)\in L^2(R)$. Problem (1.1) can now be formulated in a frequency space as follows:

$$\{\begin{array}{cc}{\xi }^2\hat{u}(\xi ,y)-{\hat{u}}_{yy}(\xi ,y)=\hat{f}(\xi ),\hfill & -\mathrm{\infty }<x<\mathrm{\infty },y>0,\hfill \\ \hat{u}(\xi ,0)=0,\hfill & -\mathrm{\infty }<x<\mathrm{\infty },\hfill \\ \hat{u}(\xi ,y){|}_{y\to \mathrm{\infty }}\text{ bounded},\hfill & -\mathrm{\infty }<x<\mathrm{\infty },\hfill \\ \hat{u}(\xi ,1)=\hat{g}(\xi ),\hfill & -\mathrm{\infty }<x<\mathrm{\infty }.\hfill \end{array}$$

(1.7)

By elementary calculations, the solution of problem (1.1) in the frequency space is given by

$$\hat{f}(\xi )=\frac{{\xi }^2}{1-e^{-\xi }}\hat{g}(\xi ),$$

(1.8)

which shows that $\hat{A}:L^2(R)\to L^2(R)$ is the multiplication operator. In addition, $\hat{A}$ is self-adjoint, i.e.,

$${\hat{A}}^{\ast }\hat{f}=\hat{A}\hat{f}=\frac{1-e^{-\xi }}{{\xi }^2}\hat{g}(\xi ).$$

(1.9)

2 Preliminaries

2.1 Dual least squares method

A general projection method for the operator equation $Af=g$, $A:X=L^2(\mathbb{R})\mapsto X=L^2(\mathbb{R})$ is generated by two subspace families $\{V_j\}$ and $\{Y_j\}$ of X, and the approximate solution $f_j\in V_j$ is defined to be the solution of the following problem:

$$〈A{f}_j,y〉=〈g,y〉,\mathrm{\forall }y\in Y_j,$$

(2.1)

where $〈\cdot ,\cdot 〉$ denotes the inner product in X. If $V_j\subset R(A^{\ast })$ and subspaces $Y_j$ are chosen such that

$$A^{\ast }{Y}_j=V_j,$$

(2.2)

then we have a special case of the projection method known as the dual least squares method. If ${\{{\psi }_{\lambda }\}}_{\lambda \in I_j}$ is an orthogonal basis of $V_j$ and $y_{\lambda }$ is the solution of the equation

$$A^{\ast }{y}_{\lambda }=k_{\lambda }{\psi }_{\lambda },\parallel y_{\lambda }\parallel =1,$$

(2.3)

then the approximate solution is explicitly given by the expression

$$f_j=\sum _{\lambda \in I_j}〈g,y_{\lambda }〉\frac{1}{k_{\lambda }}{\psi }_{\lambda }.$$

(2.4)

2.2 Subspaces $Y_j$

In this section, we investigate some properties of the subspaces $Y_j$. A method for constructing the basis of the subspace is given. This method is different from [17] in that the function $v(\xi ,y)$ is not specific. The basis of $Y_j$ cannot be explicitly obtained by dilations and integer translations of a function like the one in [17].

According to $A^{\ast }{Y}_j=V_j$, the subspaces $Y_j$ are spanned by $f_{\lambda }$, $\lambda \in I_j$, where

$$A^{\ast }{f}_{\lambda }={\mathrm{\Psi }}_{\lambda }\text{and}{k}_{\lambda }={\parallel f_{\lambda }\parallel }^{-1},y_{\lambda }=\frac{f_{\lambda }}{\parallel f_{\lambda }\parallel }=k_{\lambda }{f}_{\lambda }.$$

(2.5)

Since $supp{\hat{\mathrm{\Psi }}}_{\lambda }$ is compact, the solution exists for any $y\in [0,1)$. Similarly, the solution of the adjoint problem is unique. Therefore, for a given ${\mathrm{\Psi }}_{\lambda }$, $f_{\lambda }$ can be uniquely determined; furthermore,

$${\hat{f}}_{\lambda }=v(\xi ,y){\hat{\mathrm{\Psi }}}_{\lambda }(\xi )\iff {\hat{y}}_{\lambda }=v(\xi ,y)k_{\lambda }{\hat{\mathrm{\Psi }}}_{\lambda }(\xi ),\lambda =\{j,k\}.$$

(2.6)

As for some properties of $k_{\lambda }$, the results are similar to Lemma 3.2 of [17]. Here we omit it.

2.3 Meyer wavelets

The Meyer wavelet ψ is a function $C^{\mathrm{\infty }}(\mathbb{R})$ defined by its Fourier transform as follows [18]:

$$\hat{\psi }(\xi )=\{\begin{array}{cc}\frac{1}{\sqrt{2\pi }}{e}^{i\frac{\xi }{2}}sin[\frac{\pi }{2}\nu (\frac{3}{2\pi }|\xi |-1)],\hfill & \frac{2\pi }{3}\le |\xi |\le \frac{4\pi }{3},\hfill \\ \frac{1}{\sqrt{2\pi }}{e}^{i\frac{\xi }{2}}cos[\frac{\pi }{2}\nu (\frac{3}{4\pi }|\xi |-1)],\hfill & \frac{4\pi }{3}\le |\xi |\le \frac{8\pi }{3},\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$

(2.7)

where $\nu \in C^k$ is equal to 0 for $x\le 0$, is equal to 1 for $x\ge 1$, and $\nu (x)+\nu (1-x)=1$ for $0<x<1$. The corresponding scaling function ϕ is defined by

$$\hat{\varphi }(\xi )=\{\begin{array}{cc}\frac{1}{2\pi },\hfill & |\xi |\le \frac{2\pi }{3},\hfill \\ \frac{1}{\sqrt{2\pi }}cos[\frac{\pi }{2}\nu (\frac{3}{2\pi }|\xi |-1)],\hfill & \frac{2\pi }{3}\le |\xi |\le \frac{4\pi }{3},\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$

(2.8)

Let us list some notations: ${\psi }_{j,k}(x):=2^{\frac{j}{2}}\psi (2^{j}x-k)$, ${\varphi }_{j,k}(x):=2^{\frac{j}{2}}\varphi (2^{j}x-k)$, $j,k\in \mathbb{Z}$; ${\mathrm{\Psi }}_{-1,k}:={\varphi }_{0,k}$ and ${\mathrm{\Psi }}_{l,k}:={\psi }_{l,k}$ for $l\ge 0$; wavelet spaces $W_j=span{\{{\psi }_{j,k}\}}_{j,k\in {\mathbb{Z}}^2}$; some index sets (where $J\ge 1$ is a fixed integer)

$$\begin{array}{r}I=\{\{j,k\}:j,k\in \mathbb{Z}\}\subset {\mathbb{Z}}^2,\\ I_J=\{\{j,k\}:j=-1,0,\dots ,J-1;k\in \mathbb{Z}\}\subset {\mathbb{Z}}^2,\\ I_{j\ge J+1}=\{\{j,k\}:j\ge J;k\in \mathbb{Z}\}\subset {\mathbb{Z}}^2.\end{array}$$

(2.9)

By successively decomposing the scaling space $V_J$, $V_{J-1}$ and so on, we have $V_J=V_{J-1}\oplus W_{J-1}=V_{J-2}\oplus W_{J-2}\oplus W_{J-1}=\cdots =V_0\oplus W_1\oplus \cdots \oplus W_{J-1}$, hence we can define the subspaces $V_J$

$$V_J=\overline{span{\{{\mathrm{\Psi }}_{\lambda }\}}_{\lambda \in I_J}}.$$

(2.10)

Define an orthogonal projection $P_J:L^2(\mathbb{R})\mapsto V_J$:

$$P_J\phi =\sum _{\lambda \in I_J}〈\phi ,{\mathrm{\Psi }}_{\lambda }〉{\mathrm{\Psi }}_{\lambda },\mathrm{\forall }\phi \in L^2(\mathbb{R}),$$

(2.11)

then replace the ${\{{\psi }_{\lambda }\}}_{\lambda \in I_j}$ in (2.4) by ${\{{\mathrm{\Psi }}_{\lambda }\}}_{\lambda \in I_J}$. We easily conclude

$$f_J=P_{J}f.$$

(2.12)

From the point of view of an application to problem (1.1), the important property of Meyer wavelets is the compactness of their support in the frequency space. Indeed, since

$${\hat{\psi }}_{j,k}(\xi )=2^{-\frac{j}{2}}{e}^{-i{2}^{-j}k\xi }\hat{\psi }\left(2^{-j}\xi \right),{\hat{\varphi }}_{j,k}(\xi )=2^{-\frac{j}{2}}{e}^{-i{2}^{-j}k\xi }\hat{\varphi }\left(2^{-j}\xi \right),$$

it follows that for any $k\in \mathbb{Z}$,

$$supp({\hat{\psi }}_{j,k})=\{\xi :\frac{2}{3}\pi 2^j\le |\xi |\le \frac{8}{3}\pi 2^j\},supp({\hat{\varphi }}_{j,k})=\{\xi :|\xi |\le \frac{4}{3}\pi 2^j\}.$$

(2.13)

From (2.11), $P_J$ can be seen as a low pass filter. The frequencies with greater than $\frac{8}{3}\pi 2^J$ are filtered away.

3 Error estimates for the wavelet dual least square method

Theorem 3.1 If $f(x)$ is the solution of problem (1.1) satisfying condition (1.3), then it holds that

$$\parallel f(\cdot )-P_{J}f(\cdot )\parallel \le {\left(\frac{2}{3}\pi 2^{J+1}\right)}^{-p}E.$$

(3.1)

Proof From (2.11), we have

$$\begin{array}{r}f(\cdot )=\sum _{\lambda \in I}〈f(\cdot ),{\mathrm{\Psi }}_{\lambda }〉{\mathrm{\Psi }}_{\lambda },\\ P_{J}f(\cdot )=\sum _{\lambda \in I_J}〈f(\cdot ),{\mathrm{\Psi }}_{\lambda }〉{\mathrm{\Psi }}_{\lambda }.\end{array}$$

By virtue of the Parseval relation, with the ${\hat{\mathrm{\Psi }}}_{\lambda }$’s compact support (2.13), there holds

$$\begin{array}{rl}\parallel f(\cdot )-P_{J}f(\cdot )\parallel & =\parallel \hat{f}(\cdot )-\widehat{P_{J}f}(\cdot )\parallel =\parallel \sum _{\lambda \in I}〈\hat{f},{\hat{\mathrm{\Psi }}}_{\lambda }〉{\hat{\mathrm{\Psi }}}_{\lambda }-\sum _{\lambda \in I_J}〈\hat{f},{\hat{\mathrm{\Psi }}}_{\lambda }〉{\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ =\parallel \sum _{\lambda \in I_{j\ge J+1}}〈\hat{f},{\hat{\mathrm{\Psi }}}_{\lambda }〉{\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ =\parallel \sum _{\lambda \in I_{j\ge J+1}}〈{(1+{\xi }^2)}^{-p/2}{(1+{\xi }^2)}^{p/2}\hat{f}(\cdot ),{\hat{\mathrm{\Psi }}}_{\lambda }〉{\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ \le \underset{\frac{2}{3}\pi 2^{j+1}\le |\xi |\le \frac{8}{3}\pi 2^{j+1}}{sup}{(1+{\xi }^2)}^{-p/2}|\cdot \parallel \sum _{\lambda \in I_{j\ge J+1}}〈{(1+{\xi }^2)}^{p/2}\hat{f}(\cdot ),{\hat{\mathrm{\Psi }}}_{\lambda }〉{\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ \le {\left(\frac{2}{3}\pi 2^{J+1}\right)}^{-p}E.\end{array}$$

The approximate solution for noisy data $g_{\delta }$ is explicitly given by

$$P_J{f}^{\delta }(x)=f_J^{\delta }=\sum _{\lambda \in I_J}〈f^{\delta },{\mathrm{\Psi }}_{\lambda }〉{\mathrm{\Psi }}_{\lambda }=\sum _{\lambda \in I_J}〈g_{\delta },y_{\lambda }〉\frac{1}{k_{\lambda }}{\mathrm{\Psi }}_{\lambda }.$$

(3.2)

Now we will estimate the error $\parallel P_J{f}^{\delta }-P_{J}f\parallel$. □

Theorem 3.2 If $g_{\delta }$ is noisy data satisfying condition (1.2), then for any fixed y, we have

$$\parallel P_J{f}^{\delta }-P_{J}f\parallel \le C{\left(\frac{2}{3}\pi 2^{J+1}\right)}^2\delta ,C=\frac{4}{1-e^{-\frac{8}{3}\pi }}\doteq 4.$$

(3.3)

Proof Using (2.4), (3.6), and (4.2), from the Parseval relation, we have

$$\begin{array}{r}\parallel P_J{f}^{\delta }-P_{J}f\parallel \\ =\parallel \sum _{\lambda \in I_J}\frac{1}{1-\lambda }({\hat{g}}_{\delta }-\hat{g},\frac{{\xi }^2}{1-e^{-\xi }}{k}_{\lambda }{\hat{\mathrm{\Psi }}}_{\lambda }){\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ \le \underset{\frac{2}{3}\pi 2^J\le |\xi |\le \frac{8}{3}\pi 2^J}{sup}\frac{{\xi }^2}{1-e^{-\xi }}\parallel \sum _{\lambda \in I_J}({\hat{g}}_{\delta }-\hat{g},{\hat{\mathrm{\Psi }}}_{\lambda }){\hat{\mathrm{\Psi }}}_{\lambda }\parallel \\ \le \underset{\frac{2}{3}\pi 2^J\le |\xi |\le \frac{8}{3}\pi 2^J}{sup}\frac{{\xi }^2}{1-e^{-\xi }}\parallel \widehat{P_J(g_{\delta }-g)}\parallel \\ \le \underset{\frac{2}{3}\pi 2^J\le |\xi |\le \frac{8}{3}\pi 2^J}{sup}\frac{{\xi }^2}{1-e^{-\xi }}\parallel (g_{\delta }-g)\parallel \\ \le \frac{{(\frac{8}{3}\pi 2^J)}^2}{1-e^{-\frac{8}{3}\pi 2^J}}\delta \\ \le C{\left(\frac{2}{3}\pi 2^{J+1}\right)}^2\delta ,C=\frac{4}{1-e^{-\frac{8}{3}\pi }}\doteq 4.\end{array}$$

 □

Theorem 3.3 If $f(\cdot )$ is the solution of problem (1.1) satisfying the condition ${\parallel u(\cdot ,1)\parallel }_p\le E$, let $P_J{f}^{\delta }$ be given by (3.2). If $\parallel g-g_{\delta }\parallel \le \delta$ and $J=J(\delta )$ is selected such that

$$\frac{2}{3}\pi 2^{J+1}={\left(\frac{E}{\delta }\right)}^{\frac{1}{p+2}},$$

(3.4)

then

$$\parallel f(\cdot )-P_J{f}^{\delta }(\cdot )\parallel \le (C+1)E^{\frac{2}{p+2}}{\delta }^{\frac{p}{p+2}}.$$

(3.5)

Proof By the triangle inequality, we have

$$\parallel f(\cdot )-P_J{f}^{\delta }(\cdot )\parallel \le \parallel f(\cdot )-P_{J}f(\cdot )\parallel +\parallel P_{J}f(\cdot )-P_J{f}^{\delta }(\cdot )\parallel .$$

(3.6)

Combining Theorem 3.1 with Theorem 3.2, we obtain the convergence estimate of our method. □

4 Numerical tests

Example 1 It is easy to see that the function

$$u(x,y)=(1-e^{-y})sinx$$

(4.1)

and the function

$$f(x)=sinx$$

(4.2)

satisfy problem (1.1) with exact data

$$g(x)=(1-e^{-1})sinx.$$

(4.3)

We will do the numerical tests in the interval $x\in [-10,10]$.

Suppose that the sequence $g{(x_i)}_{i=1}^n$ represents samples from the function $g(x)$ on an equidistant grid, and n is even, then we add a random uniformly distributed perturbation to each data, and obtain the perturbation data

$$g_{\delta }=g+\mu \ast rand(size(g)),$$

(4.4)

where

$$g=(g(x_1),g(x_2),\dots ,g(x_n)),x_i=(i-1)\mathrm{\Delta }x-10,\mathrm{\Delta }x=\frac{20}{n-1},i=1,2,\dots ,n.$$

Then the total noise δ can be measured in the sense of root mean square error according to

$$\delta ={\parallel g_{\delta }-g\parallel }_{l^2}={\left(\frac{1}{n}\sum _{i=1}^n{({(g_{\delta })}_i-g_i)}^2\right)}^{1/2}.$$

(4.5)

The computed errors are defined by

$$l^2\text{ norm error: }E(f)=\sqrt{\frac{1}{n}\sum _{i=1}^n{(f(x_i)-f_r(x_i))}^2},$$

(4.6)

$$\text{relative error: }ER(f)=\sqrt{\frac{1}{n}\sum _{i=1}^n{(f(x_i)-f_r(x_i))}^2}/\sqrt{\frac{1}{n}\sum _{i=1}^n(f{(x_i)}^2},$$

(4.7)

where $x_i$, $i=1,2,\dots ,n$, are the test points. In computation, we take $n=64$, $f_r(\cdot )$ denotes the regularization solution. In numerical tests, the regularization parameter α is selected by $\alpha =\frac{{\delta }^2}{M^2}$.

From Figures 1 and 2, we can conclude that the approximation effect of wavelet dual least square regularization for $y=1$ and $y=0.5$.

Figure 1

Data φ .

Figure 2

Noisy data ${\mathit{\phi }}_{\mathit{\delta }}$ .