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

Authors: Xiao-Xiao Li; Fan Yang; Jie Liu; Lan Wang
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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(u(x,y),f(x)) which satisfy(1)▵u(x,y)-k2u(x,y)=f(x),-∞<x<∞,0<y<+∞,u(x,0)=0,-∞<x<∞,u(x,y)|y→∞bounded,-∞<x<∞,u(x,1)=g(x),-∞<x<∞,
where f(x) is the unknown source depending only on one spatial variable, u(x,1)=g(x) is the supplementary condition, and the constant k>0 is the wave number. In applications, input data g(x) can only be measured; there will be measured data function gδ(x) which is merely in L2(ℝ) and satisfies
(2)∥g-gδ∥L2(ℝ)≤δ,
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 f^(ξ) denote the Fourier transform of f(x)∈L2(ℝ) which is defined by
(3)f^(ξ):=12π∫-∞∞‍e-iξxf(x)dx.
The problem (1) can now be formulated in frequency space as follows:
(4)u^yy(ξ,y)-(ξ2+k2)u^(ξ,y)=f^(ξ),ξ∈ℝ,y>0,u^(ξ,0)=0,ξ∈ℝ,u^(ξ,y)|y→∞bounded,ξ∈ℝ,u^(ξ,1)=g^(ξ),ξ∈ℝ.
The solution of the problem (4) is given by
(5)f^(ξ)=-ξ2+k21-e-ξ2+k2g^(ξ).
So,
(6)f(x)=-12π∫-∞∞‍eiξxξ2+k21-e-ξ2+k2g^(ξ)dξ.
The unbounded function (ξ2+k2)/(1-e-ξ2+k2) in (5) or (6) can be seen as an amplification factor of g^(ξ) when ξ→∞. Therefore, when we consider our problem in L2(ℝ), the exact data function g^(ξ) must decay. But in the applications, the input data g(x) can only be measured and never be exact. Thus, if we try to obtain the unknown source f(x), 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,
(7)∥f(·)∥Hp≤E,p>0,
where E>0 is a constant; ∥·∥Hp denotes the norm in Sobolev space Hp(ℝ) defined by
(8)∥f(·)∥Hp:=(∫-∞∞‍|f^(ξ)|2(1+ξ2)pdξ)1/2.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)▵u(x,y)-k2u(x,y)+μ2fxx(x)=f(x),
where the choice of μ is based on some a priori knowledge about the magnitude of the errors in the data gδ. 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. Section2 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.
Forr≥1, there holds
(10)11-e-r<2.Lemma 2.
For0<μ<1, there hold the following inequalities:
(11)supξ∈ℝ(1-11+ξ2μ2)(1+ξ2)-(p/2)≤max{μp,μ2},(12)supξ∈ℝξ2+k2(1+ξ2μ2)(1-e-ξ2+k2)≤2μ2+2k2.Proof.
Let(13)G(ξ):=(1-11+ξ2μ2)(1+ξ2)-(p/2).
The proof of (11) can be separated from three cases.
Case  1
(
|
ξ
|
≥
ξ
0
:
=
1
/
μ
). We get
(14)G(ξ)≤(1+ξ2)-(p/2)≤|ξ|-p≤ξ0-p=μp.
Case  2
(
1
<
|
ξ
|
<
ξ
0
). We obtain
(15)G(ξ)=ξ2μ21+ξ2μ2(1+ξ2)-(p/2)≤ξ2-pμ21+ξ2μ2≤ξ2-pμ2.
If 0<p≤2, we have
(16)G(ξ)<ξ02-pμ2=μp.
If p>2, then
(17)G(ξ)<μ2.
Case  3
(
|
ξ
|
≤
1
). We get
(18)G(ξ)≤ξ2μ2(1+ξ2)-(p/2)≤ξ2μ2≤μ2.
Now combining (14) with (16), (17), and (18), the first inequality (11) holds.
Let(19)B(ξ):=ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2),D(ξ):=ξ2+k21-e-ξ2+k2.
Like the above proof, the proof of (12) is divided into two cases.
Case  1
(
|
ξ
|
≤
ξ
0
:
=
1
/
μ
). Using Lemma 1, we have
(20)D(ξ)≤D(1μ)≤2μ2+2k2.
So,
(21)B(ξ)≤D(ξ)≤2μ2+2k2.
Case  2
(
|
ξ
|
>
ξ
0
)
. Using Lemma 1, we get
(22)B(ξ)=ξ2(1-e-ξ2+k2)(1+ξ2μ2)+k2(1-e-ξ2+k2)(1+ξ2μ2)(23)≤2ξ21+ξ2μ2+k21-e-ξ2+k2≤2μ2+2k2.
Combining (21) with (23), the second inequality (12) holds.


## 3. A Quasireversibility Regularization Method

In this section, we consider the following system:(24)▵u(x,y)-k2u(x,y)+μ2fxx(x)=f(x),Δuu(x,y)-∞<x<∞,0<y<+∞,u(x,0)=0,-∞<x<∞,u(x,y)|y→∞bounded,-∞<x<∞,u(x,1)=gδ(x),-∞<x<∞,
where the parameter μ is regarded as a regularization parameter. The problem (24) can be formulated in frequency space as follows:
(25)u^yy(ξ,y)-(ξ2+k2)u^(ξ,y)=(1+ξ2μ2)f^(ξ),u^yy(ξ,y)-(ξ2+k2)u^(ξ,y)=ξ∈ℝ,y>0,u^(ξ,0)=0,ξ∈ℝ,u^(ξ,y)|y→∞bounded,ξ∈ℝ,u^(ξ,1)=g^δ(ξ),ξ∈ℝ.
The solution to the problem (25) is given by
(26)f^(ξ)=-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^δ(ξ):=f^μδ(ξ),
so
(27)fμδ(x)=-12π∫-∞∞‍eiξxξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^δ(ξ)dξ,
which is called the quasireversibility regularization solution.It is easy to see that, for smallμ, when |ξ| is small, (ξ2+k2)/(1+ξ2μ2) is close to ξ2+k2. When |ξ| becomes large, (ξ2+k2)/(1+ξ2μ2) is bounded. So, fμδ(x) is considered as an approximation of f(x).Now we will give a convergence error estimate between the regularization solution and the exact solution by the following theorem.Theorem 3.
Letf(x) be the solution of (1) whose Fourier transform is given by (5). Let gδ(x) be the measured data at y=1 satisfying (2).  Let priori condition (7) hold for p>0. Let fμδ(x) be the quasireversibility regularization approximation to f(x) given by (27). If one selects
(28)μ=(δE)1/(p+2),
then one obtains the following error estimate:
(29)∥f(·)-fμδ(·)∥≤2δp/(p+2)E2/(p+2)(1+12max{1,(δE)(2-p)/(p+2)})+2k2δ.Proof.
Due to the Parseval formula and the triangle inequality, we have(30)∥f(·)-fμδ(·)∥=∥f^(·)-f^μδ(·)∥=∥(-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)-ξ2+k21-e-ξ2+k2g^(ξ)-(-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^δ(ξ))∥=∥ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)ξ2+k21-e-ξ2+k2g^(ξ)-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^δ(ξ)∥≤∥ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)ξ2+k21-e-ξ2+k2g^(ξ)-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^(ξ)∥+∥ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^(ξ)-ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)g^δ(ξ)∥=∥ξ2+k21-e-ξ2+k2g^(ξ)(1-11+ξ2μ2)∥+∥ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2)(g^(ξ)-g^δ(ξ))∥≤∥f^(ξ)(1+ξ2)p/2(1-11+ξ2μ2)(1+ξ2)-(p/2)∥+supξ∈ℝ(ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2))×∥g^(ξ)-g^δ(ξ)∥≤supξ∈ℝ((1-11+ξ2μ2)(1+ξ2)-(p/2))×∥f^(ξ)(1+ξ2)p/2∥+supξ∈ℝ(ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2))δ≤supξ∈ℝ((1-11+ξ2μ2)(1+ξ2)-(p/2))E+supξ∈ℝ(ξ2+k2(1-e-ξ2+k2)(1+ξ2μ2))δ≤max{μp,μ2}E+2μ2δ+2k2δ=max{(δE)p/(p+2),(δE)2/(p+2)}E+2(δE)-2/(p+2)δ+2k2δ=2δp/(p+2)E2/(p+2)(1+12max{1,(δE)(2-p)/(p+2)})+2k2δ.Remark 4.
If0<p≤2,
(31)∥f(·)-fμδ(·)∥≤3δp/(p+2)E2/(p+2)+2k2δ→0asδ→0.
If p>2,
(32)∥f(·)-fμδ(·)∥≤2δp/(p+2)E2/(p+2)+δ2/(p+2)Ep/(p+2)+2k2δ→0asδ→0.
Hence, fμδ(x) can be viewed as the approximation of the exact solution f(x).


## 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 solutionf(x) and obtained the exact data function g(x) through solving the forward problems. Then we added a normally distributed perturbation to each data function and obtained vectors gδ(x). 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(u(x,y),f(x)).Example 5.
It is easy to see that the function(33)u(x,y)={(1-e-2ky)sinkx,y>0,0,y≤0
and the function
(34)f(x)=-2k2sinkx
are satisfied with the problem (1) with exact data
(35)g(x)=(1-e-2k)sinkx.
Suppose that the sequence g(xi)i=1i=n represents samples from the function g(x) on an equidistant grid, then we use the rand function given in MATLAB to generate the noisy data,
(36)(gδ)=g+εrand(size(g)),
where
(37)g=(g(x1),…,g(xn))T,xi=(i-1)Δx-10,Δx=20n-1,i=1,…,n.
The function “rand(·)” generates arrays of random numbers whose elements are normally distributed with mean 0, variance σ2=1. “Rand(size(g))” returns an array of random entries that is the same size as g. The total noise level δ can be measured in the sense of Root Mean Square Error (RMSE) according to
(38)δ=∥gδ-g∥2=(1n∑i=1n‍(gi-giδ)2)1/2.In our computations, we taken=100, and the relative error is given as follows:
(39)rerr(f):=∥fμδ-f∥2∥f∥2,
where ∥·∥2 is defined by (38).Tables1, 2, 3 and 4 show that parameters δ, μ, and rerr(f) all depend on the perturbation ε. Parameters δ, μ and rerr(f) decrease with the decrease of ε. These are consistent with our error estimate. In addition, rerr(f) does not decrease for stronger “smoothness” assumptions on the exact solution f(x).Table 1
δ, μ and the relative error rerr(f) with k=2, p=1 and μ=(δ/E)1/3 for different ε.
ε
1
0
-
1
1
0
-
2
1
0
-
3
1
0
-
4
1
0
-
5
δ
0.1866
0.0205
0.0020
2.0115
×
1
0
-
4
1.9050
×
1
0
-
5
μ
0.0727
0.0389
0.0177
0.0083
0.0038
rerr(f)
2.0442
1.7760
0.7874
0.5816
0.4817Table 2
δ, μ and the relative error rerr(f) with k=2, p=2 and μ=(δ/E)1/4 for different ε.
ε
1
0
-
1
1
0
-
2
1
0
-
3
1
0
-
4
1
0
-
5
δ
0.2119
0.0192
0.0020
1.9216
×
1
0
-
4
1.994
×
1
0
-
5
μ
0.0080
0.0555
0.0318
0.0177
0.0100
rerr(f)
1.9593
1.6765
0.6422
0.4318
0.3215Table 3
δ, μ and the relative error rerr(f) with k=2, p=3 and μ=(δ/E)1/5 for different ε.
ε
1
0
-
1
1
0
-
2
1
0
-
3
1
0
-
4
1
0
-
5
δ
0.2026
0.0207
0.0020
2.0160
×
1
0
-
4
1.9794
×
1
0
-
5
μ
0.0813
0.0615
0.0387
0.0244
0.0153
rerr(f)
1.7544
1.6173
0.5137
0.3176
0.2689Table 4
δ, μ and the relative error rerr(f) with k=2, p=5 and μ=(δ/E)1/7 for different ε.
ε
1
0
-
1
1
0
-
2
1
0
-
3
1
0
-
4
1
0
-
5
δ
0.2085
0.0206
0.0021
2.0920
×
1
0
-
4
2.0210
×
1
0
-
5
μ
0.0770
0.0639
0.0464
0.0334
0.0240
rerr(f)
2.0533
1.6183
0.6187
0.4769
0.3231Example 6.
Consider a piecewise smooth unknown source:(40)f(x)={0,0≤-10≤-5,x+5,-5≤x≤0,-x+5,0≤x≤5,0,5≤x≤10.Example 7.
Consider the following discontinuous unknown source:(41)f(x)={-1,0≤-10≤-5,1,-5≤x≤0,-1,0≤x≤5,1,5≤x≤10.From Figures1, 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   k 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 f(x) does not have an analytical solution, the data g(x) 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 withk=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 withk=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 withk=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 withk=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 withk=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 withk=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 withk=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 withk=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 withk=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.



---

*Source: 245963-2013-04-24.xml*