The Quasireversibility Regularization Method for Identifying the Unknown Source for the Modified Helmholtz Equation

Abstract

This paper discusses the problem of determining an unknown source which depends only on one variable for the modified Helmholtz equation. This problem is ill-posed in the sense that the solution (if it exists) does not depend continuously on the data. The regularization solution is obtained by the quasireversibility regularization method. Convergence estimate is presented between the exact solution and the regularization solution. Moreover, numerical results are presented to illustrate the accuracy and efficiency of this method.

1. Introduction

Inverse source problems arise in many branches of science and engineering, for example, heat conduction, crack identification, electromagnetic theory, geophysical prospecting, and pollutant detection. For the heat source identification, there have been a large number of research results for different forms of heat source [1–5]. The modified Helmholtz equation or the Yukawa equation which is pointed out in [6] appears in implicit marching schemes for the heat equation, in Debye-Hückel theory, and in the linearization of the Poisson-Boltzmann equation. The underlying free-space Green’s function is usually referred to as the Yukawa potential in nuclear physics. In physics, chemistry, and biology, when Coulomb forces are damped by screening effects, this Green’s function is also known as the screened Coulomb potential. To the authors’ knowledge, there were few papers for identifying the unknown source on the modified Helmholtz equation by regularization method.

In this paper, we consider the following inverse problem: to find a pair of functions (𝑢(𝑥,𝑦),𝑓(𝑥)) which satisfy▵𝑢𝑥,𝑦−𝑘²𝑢(𝑥),𝑢()𝑢𝑥,𝑦=𝑓−∞<𝑥<∞,0<𝑦<+∞,𝑥,0=0,−∞<𝑥<∞,𝑥,𝑦_𝑦→∞𝑢()(𝑥)bounded,−∞<𝑥<∞,𝑥,1=𝑔,−∞<𝑥<∞, (1) where 𝑓(𝑥) is the unknown source depending only on one spatial variable, 𝑢(𝑥,1)=𝑔(𝑥) is the supplementary condition, and the constant 𝑘>0 is the wave number. In applications, input data 𝑔(𝑥) can only be measured; there will be measured data function 𝑔^𝛿(𝑥) which is merely in 𝐿²(ℝ) and satisfies 𝑔−𝑔^𝛿_𝐿2(ℝ)≤𝛿, (2) where the constant 𝛿>0 represents a noise level of input data.

The ill-posedness can be seen by solving the problem (1) in the Fourier domain. Let 𝑓(𝜉) denote the Fourier transform of 𝑓(𝑥)∈𝐿²(ℝ) which is defined by 𝑓(𝜉)1∶=√2𝜋^∞−∞𝑒^{−𝑖𝜉𝑥}𝑓(𝑥)𝑑𝑥. (3) The problem (1) can now be formulated in frequency space as follows: 𝑢_𝑦𝑦𝜉,𝑦−𝜉²+𝑘²𝑢=𝑓(𝜉)𝑢()𝑢𝜉,𝑦,𝜉∈ℝ,𝑦>0,𝜉,0=0,𝜉∈ℝ,𝜉,𝑦_𝑦→∞𝑢()=𝑔(𝜉)bounded,𝜉∈ℝ,𝜉,1,𝜉∈ℝ. (4) The solution of the problem (4) is given by 𝑓(𝜉)𝜉=−²+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉). (5) So, 𝑓(𝑥)1=−√2𝜋^∞−∞𝑒^𝑖𝜉𝑥𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉)𝑑𝜉. (6) The unbounded function (𝜉²+𝑘²)/(1−𝑒^{−𝜉2+𝑘2}) in (5) or (6) can be seen as an amplification factor of 𝑔(𝜉) when 𝜉→∞. Therefore, when we consider our problem in 𝐿²(ℝ), the exact data function 𝑔(𝜉) must decay. But in the applications, the input data 𝑔(𝑥) can only be measured and never be exact. Thus, if we try to obtain the unknown source 𝑓(𝑥), high frequency components in the error are magnified and can destroy the solution. So in the following section, we will use the regularization method to deal with the ill-posed problem. Before doing that, we impose an a priori bound on the input data, that is, 𝑓(⋅)_𝐻𝑝≤𝐸,𝑝>0, (7) where 𝐸>0 is a constant; ‖⋅‖_𝐻^𝑝 denotes the norm in Sobolev space 𝐻^𝑝(ℝ) defined by 𝑓(⋅)_𝐻𝑝∶=^∞_−∞𝑓(𝜉)²1+𝜉^2𝑝𝑑𝜉^1/2. (8)

In this paper, a new regularization method which is proposed as an alternative way of regularization methods for identifying unknown source, is given. Actually, we discuss the possibility of modifying (1) to obtain a stable approximation; that is, we will investigate the following problem: ▵𝑢𝑥,𝑦−𝑘²𝑢𝑥,𝑦+𝜇²𝑓_𝑥𝑥(𝑥)=𝑓(𝑥), (9) where the choice of 𝜇 is based on some a priori knowledge about the magnitude of the errors in the data 𝑔^𝛿. The idea is called the quasireversibility regularization method. We were inspired from Eldén [7] who considered a standard inverse heat conduction problem and the idea initially came from Weber [8]. Now the quasireversibility regularization method has been studied for solving various types of inverse problems [8–13].

In nature, the quasireversibility regularization method transfers an ill-posed problem to an approximate well-posed problem which can be discretized using standard technique, for example, finite differences. For the numerical implementation of the quasireversibility regularization method, one can refer to [9–13]. Our aim here is to discuss the stability and convergence analysis of regularization method.

This paper is organized as follows. Section 2 gives some auxiliary results. Section 3 gives a regularization solution and error estimation. Section 4 gives three examples to illustrate the accuracy and efficiency of this method. Section 5 puts an end to this paper with a brief conclusion.

2. Some Auxiliary Results

Now we give some important lemmas, which are very useful for our main conclusion.

Lemma 1. For 𝑟≥1, there holds 11−𝑒^−𝑟<2. (10)

Lemma 2. For 0<𝜇<1, there hold the following inequalities: sup_𝜉∈ℝ11−1+𝜉²𝜇²1+𝜉^2−(𝑝/2)≤max𝜇^𝑝,𝜇²,sup_𝜉∈ℝ𝜉²+𝑘²1+𝜉²𝜇²1−𝑒^{−𝜉2+𝑘2}≤2𝜇²+2𝑘². (11)(12)

Proof. Let 𝐺(𝜉)1∶=1−1+𝜉²𝜇²1+𝜉^2−(𝑝/2). (13) The proof of (11) can be separated from three cases.
Case  1 (|𝜉|≥𝜉₀∶=1/𝜇). We get 𝐺(𝜉)≤1+𝜉^2−(𝑝/2)≤𝜉^−𝑝≤𝜉₀^−𝑝=𝜇^𝑝. (14)
Case  2 (1<|𝜉|<𝜉₀). We obtain 𝐺(𝜉)=𝜉²𝜇²1+𝜉²𝜇²1+𝜉^2−(𝑝/2)≤𝜉^2−𝑝𝜇²1+𝜉²𝜇²≤𝜉^2−𝑝𝜇². (15) If 0<𝑝≤2, we have 𝐺(𝜉)<𝜉₀^2−𝑝𝜇²=𝜇^𝑝. (16) If 𝑝>2, then 𝐺(𝜉)<𝜇². (17)
Case  3 (|𝜉|≤1). We get 𝐺(𝜉)≤𝜉²𝜇²1+𝜉^2−(𝑝/2)≤𝜉²𝜇²≤𝜇². (18) Now combining (14) with (16), (17), and (18), the first inequality (11) holds.
Let 𝐵(𝜉)𝜉∶=²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²,𝐷(𝜉)𝜉∶=²+𝑘²1−𝑒^{−𝜉2+𝑘2}. (19) Like the above proof, the proof of (12) is divided into two cases.
Case  1 (|𝜉|≤𝜉₀∶=1/𝜇). Using Lemma 1, we have 𝐷(𝜉)1≤𝐷𝜇≤2𝜇²+2𝑘². (20) So, 𝐵(𝜉)(𝜉)≤2≤𝐷𝜇²+2𝑘². (21)
Case  2 (|𝜉|>𝜉₀). Using Lemma 1, we get 𝐵(𝜉)=𝜉²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²≤2𝜉²1+𝜉²𝜇²+𝑘²1−𝑒^{−𝜉2+𝑘2}≤2𝜇²+2𝑘². (22)(23) Combining (21) with (23), the second inequality (12) holds.

3. A Quasireversibility Regularization Method

In this section, we consider the following system: ▵𝑢𝑥,𝑦−𝑘²𝑢𝑥,𝑦+𝜇²𝑓_𝑥𝑥(𝑥)(𝑥),𝑢()𝑢=𝑓−∞<𝑥<∞,0<𝑦<+∞,𝑥,0=0,−∞<𝑥<∞,𝑥,𝑦_𝑦→∞𝑢()bounded,−∞<𝑥<∞,𝑥,1=𝑔^𝛿(𝑥),−∞<𝑥<∞, (24) where the parameter 𝜇 is regarded as a regularization parameter. The problem (24) can be formulated in frequency space as follows: 𝑢_𝑦𝑦𝜉,𝑦−𝜉²+𝑘²𝑢𝜉,𝑦=1+𝜉²𝜇²𝑓(𝜉),𝑢()𝑢𝜉∈ℝ,𝑦>0,𝜉,0=0,𝜉∈ℝ,𝜉,𝑦_𝑦→∞𝑢()=𝑔bounded,𝜉∈ℝ,𝜉,1^𝛿(𝜉),𝜉∈ℝ. (25) The solution to the problem (25) is given by 𝑓(𝜉)𝜉=−²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔^𝛿(𝜉)∶=𝑓^𝛿_𝜇(𝜉), (26) so 𝑓^𝛿_𝜇(𝑥)1=−√2𝜋^∞−∞𝑒^𝑖𝜉𝑥𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔^𝛿(𝜉)𝑑𝜉, (27) which is called the quasireversibility regularization solution.

It is easy to see that, for small 𝜇, when |𝜉| is small, (𝜉²+𝑘²)/(1+𝜉²𝜇²) is close to 𝜉²+𝑘². When |𝜉| becomes large, (𝜉²+𝑘²)/(1+𝜉²𝜇²) is bounded. So, 𝑓^𝛿_𝜇(𝑥) is considered as an approximation of 𝑓(𝑥).

Now we will give a convergence error estimate between the regularization solution and the exact solution by the following theorem.

Theorem 3. Let 𝑓(𝑥) be the solution of (1) whose Fourier transform is given by (5). Let 𝑔^𝛿(𝑥) be the measured data at 𝑦=1 satisfying (2). Let priori condition (7) hold for 𝑝>0. Let 𝑓^𝛿_𝜇(𝑥) be the quasireversibility regularization approximation to 𝑓(𝑥) given by (27). If one selects 𝛿𝜇=𝐸^1/(𝑝+2), (28) then one obtains the following error estimate: 𝑓(⋅)−𝑓^𝛿_𝜇(⋅)≤2𝛿^{𝑝/(𝑝+2)}𝐸^2/(𝑝+2)11+2𝛿max1,E^{(2−𝑝)/(𝑝+2)}+2𝑘²𝛿. (29)

Proof. Due to the Parseval formula and the triangle inequality, we have 𝑓(⋅)−𝑓^𝛿_𝜇(⋅=𝑓(⋅)−𝑓)_𝛿𝜇(⋅=−𝜉₂)+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉)−−𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔^𝛿(𝜉)=𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉)−𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔^𝛿(𝜉)≤𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉)−𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔(𝜉)+𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔(𝜉)−𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔^𝛿(𝜉)=𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}𝑔(𝜉1)1−1+𝜉²𝜇²+𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝑔(𝜉)−𝑔^𝛿(𝜉)≤𝑓(𝜉)1+𝜉^2𝑝/211−1+𝜉²𝜇²1+𝜉^2−(𝑝/2)+sup_𝜉∈ℝ𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²×𝑔(𝜉)−𝑔^𝛿(𝜉)≤sup_𝜉∈ℝ11−1+𝜉²𝜇²1+𝜉^2−(𝑝/2)×𝑓(𝜉)1+𝜉^2𝑝/2+sup_𝜉∈ℝ𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝛿≤sup_𝜉∈ℝ11−1+𝜉²𝜇²1+𝜉^2−(𝑝/2)𝐸+sup_𝜉∈ℝ𝜉²+𝑘²1−𝑒^{−𝜉2+𝑘2}1+𝜉²𝜇²𝛿≤max𝜇^𝑝,𝜇²2𝐸+𝜇²𝛿+2𝑘²𝛿𝛿=max𝐸^{𝑝/(𝑝+2)},𝛿𝐸^2/(𝑝+2)𝐸𝛿+2𝐸^{−2/(𝑝+2)}𝛿+2𝑘²𝛿=2𝛿^{𝑝/(𝑝+2)}𝐸^2/(𝑝+2)11+2𝛿max1,𝐸^{(2−𝑝)/(𝑝+2)}+2𝑘²𝛿. (30)

Remark 4. If 0<𝑝≤2, 𝑓(⋅)−𝑓^𝛿_𝜇(⋅)≤3𝛿^{𝑝/(𝑝+2)}𝐸^2/(𝑝+2)+2𝑘²𝛿→0as𝛿→0. (31) If 𝑝>2, 𝑓(⋅)−𝑓^𝛿_𝜇(⋅)≤2𝛿^{𝑝/(𝑝+2)}𝐸^2/(𝑝+2)+𝛿^2/(𝑝+2)𝐸^{𝑝/(𝑝+2)}+2𝑘²𝛿→0as𝛿→0. (32) Hence, 𝑓^𝛿_𝜇(𝑥) can be viewed as the approximation of the exact solution 𝑓(𝑥).

4. Numerical Example

In this section, We will give three different type examples to verify the validity of the theoretical results of this method.

The numerical examples were constructed in the following way: First we selected the exact solution 𝑓(𝑥) and obtained the exact data function 𝑔(𝑥) through solving the forward problems. Then we added a normally distributed perturbation to each data function and obtained vectors 𝑔^𝛿(𝑥). Finally we obtained the regularization solutions through solving the inverse problem.

In the following, we first give an example which has the exact expression of the solutions (𝑢(𝑥,𝑦),𝑓(𝑥)).

Example 5. It is easy to see that the function 𝑢𝑥,𝑦=1−𝑒^{−√2𝑘𝑦}sin𝑘𝑥,𝑦>0,0,𝑦≤0 (33) and the function 𝑓(𝑥)=−2𝑘²sin𝑘𝑥 (34) are satisfied with the problem (1) with exact data 𝑔(𝑥)=1−𝑒^−√2𝑘sin𝑘𝑥. (35) Suppose that the sequence 𝑔(𝑥_𝑖)^𝑖=𝑛_𝑖=1 represents samples from the function 𝑔(𝑥) on an equidistant grid, then we use the rand function given in MATLAB to generate the noisy data, 𝑔^𝛿=𝑔+𝜀randsize𝑔, (36) where 𝑔=𝑔𝑥₁,…,𝑔𝑥_𝑛^𝑇,𝑥_𝑖=(𝑖−1)Δ𝑥−10,Δ𝑥=20𝑛−1,𝑖=1,…,𝑛. (37) The function “rand(⋅)” generates arrays of random numbers whose elements are normally distributed with mean 0, variance 𝜎²=1. “Rand(size(𝑔))” returns an array of random entries that is the same size as 𝑔. The total noise level 𝛿 can be measured in the sense of Root Mean Square Error (RMSE) according to 𝛿=𝑔^𝛿−𝑔₂=1𝑛^𝑛_𝑖=1𝑔_𝑖−𝑔^𝛿_𝑖^21/2. (38)

In our computations, we take 𝑛=100, and the relative error is given as follows: 𝑓𝑓rerr∶=^𝛿𝜇−𝑓₂𝑓₂, (39) where ‖⋅‖₂ is defined by (38).

Tables 1, 2, 3 and 4 show that parameters 𝛿, 𝜇, and rerr(𝑓) all depend on the perturbation 𝜀. Parameters 𝛿, 𝜇 and rerr(𝑓) decrease with the decrease of 𝜀. These are consistent with our error estimate. In addition, rerr(𝑓) does not decrease for stronger “smoothness” assumptions on the exact solution 𝑓(𝑥).

Table 1 

𝛿, 𝜇 and the relative error rerr(𝑓) with 𝑘=2, 𝑝=1 and 𝜇=(𝛿/𝐸)^1/3 for different 𝜀.

Table 2 

𝛿, 𝜇 and the relative error rerr(𝑓) with 𝑘=2, 𝑝=2 and 𝜇=(𝛿/𝐸)^1/4 for different 𝜀.

Table 3 

𝛿, 𝜇 and the relative error rerr(𝑓) with 𝑘=2, 𝑝=3 and 𝜇=(𝛿/𝐸)^1/5 for different 𝜀.

Table 4 

𝛿, 𝜇 and the relative error rerr(𝑓) with 𝑘=2, 𝑝=5 and 𝜇=(𝛿/𝐸)^1/7 for different 𝜀.

Example 6. Consider a piecewise smooth unknown source: 𝑓(𝑥)=0,0≤−10≤−5,𝑥+5,−5≤𝑥≤0,−𝑥+5,0≤𝑥≤5,0,5≤𝑥≤10. (40)

Example 7. Consider the following discontinuous unknown source: 𝑓(𝑥)=−1,0≤−10≤−5,1,−5≤𝑥≤0,−1,0≤𝑥≤5,1,5≤𝑥≤10. (41)

From Figures 1, 2, 3, 4, 5, 6, 7, 8, and 9, we can find that the smaller 𝜀 is, the better the computed approximation is and the smaller 𝑘 is, the better the computed approximation is. This is consistent with (29). From Figures 4–9, it can be seen that the numerical solution is less ideal than that of Example 5. In Examples 6 and 7, since the direct problem with the source term 𝑓(𝑥) does not have an analytical solution, the data 𝑔(𝑥) is obtained by solving the direct problem. It is not difficult to see that the well-known Gibbs phenomenon and the recovered data near the nonsmooth and discontinuities points are not accurate. Taking into consideration of the ill-posedness of the problem, the results presented in Figures 4–9 are reasonable.

Figure 1 

Comparison between the exact solution and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.1, 𝜀=0.01, 𝜀=0.001 for Example 5.

Figure 2 

Comparison between the exact solution and its computed approximations with 𝑘=2 and various noise levels of 𝜀=0.01, 𝜀=0.001, 𝜀=0.0001 for Example 5.

Figure 3 

Comparison between the exact solution (- - -) and its computed approximations with 𝑘=3 and various noise levels of 𝜀=0.01(−∗−), 𝜀=0.001(−⊳−), 𝜀=0.0001(−o−) for Example 5.

Figure 4 

Comparison between the exact solution and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.1, 𝜀=0.01, 𝜀=0.001 for Example 6.

Figure 5 

Comparison between the exact solution and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.1, 𝜀=0.01, 𝜀=0.001 for Example 6.

Figure 6 

Comparison between the exact solution and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.1, 𝜀=0.01, 𝜀=0.001 for Example 6.

Figure 7 

Comparison between the exact solution (- - -) and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.01(−∗−), 𝜀=0.001(−⊳−), 𝜀=0.0001(−o−) for Example 7.

Figure 8 

Comparison between the exact solution (- - -) and its computed approximations with 𝑘=2 and various noise levels of 𝜀=0.01(−∗−), 𝜀=0.001(−⊳−), 𝜀=0.0001(−o−) for Example 7.

Figure 9 

Comparison between the exact solution (- - -) and its computed approximations with 𝑘=1 and various noise levels of 𝜀=0.01(−∗−), 𝜀=0.001(−⊳−), 𝜀=0.0001(−o−) for Example 7.

5. Conclusions

In this paper, we considered the inverse problem of determining the unknown source using the quasireversibility regularization method for the modified Helmholtz equation. It was shown that, with a certain choice of the parameter, a stability estimate was obtained. Meanwhile, the numerical example verified the efficiency and accuracy of this method.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and helpful suggestions to improve the earlier version of the paper. The project is supported by the National Natural Science Foundation of China (no. 11171136 and no. 11261032), the Distinguished Young Scholars Fund of Lanzhou University of Technology (Q201015), and the basic scientific research business expenses of Gansu Province College.