# An Iterative Regularization Method to Solve the Cauchy Problem for the Helmholtz Equation

Authors: Hao Cheng; Ping Zhu; Jie Gao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A regularization method for solving the Cauchy problem of the Helmholtz equation is proposed. Thea priori and a posteriori rules for choosing regularization parameters with corresponding error estimates between the exact solution and its approximation are also given. The numerical example shows the effectiveness of this method.


---

## Body

## 1. Introduction

The Cauchy problem for the Helmholtz equation arises naturally in many areas of engineering and science, especially in wave propagation and vibration phenomena, such as the vibration of a structure [1], the acoustic cavity problem [2], the radiation wave [3], and the scattering of a wave [4]. However, this problem is severely ill-posed in the sense that a small change in the Cauchy data would lead to a dramatic variation in the solution. Therefore, it is necessary to study different highly efficient algorithms to solve this problem. Recently, a few special numerical methods to deal with this problem have been developed, such as the boundary element method [5], the method of fundamental solutions [6], the conjugate gradient method [7], the Landweber method [8], wavelet moment method [9], quasi-reversibility and truncation methods [10], modified Tikhonov regularization method [11, 12], the fourier regularization method [13], and so forth [9, 14, 15]. However, most of them choose the regularization parameter by thea priori rule, which depends seriously on thea priori bound E. However, in general, thea priori bound E cannot be known exactly in practice, and working with a wrong constant E may lead to a bad regularized solution. Therefore, giving thea posteriori parameter choice rule is a very meaningful topic.In this paper we will consider the following problem with inhomogeneous Dirichlet data in a strip domain:(1)Δu(x,y)+k2u(x,y)=0,x∈(0,1),y∈ℝ,u(0,y)=φ(y),y∈ℝ,ux(0,y)=0,y∈ℝ,
where the constant k>0 is the number of wave. The solution u(x,y) for 0<x<1 will be determined from the noisy data φδ(y). In this paper a regularization method of iteration type for solving this problem will be given. By dint of this method, thea priori anda posteriori rule for choosing a regularization parameter with strict theory analysis, as well as order optimal error estimates, will be obtained.The outline of the paper is as follows. In Section2, an order optimal error estimate is obtained for thea priori parameter choice rule. Thea posteriori parameter choice rule is given in Section 3, which also leads to a Hölder-type error estimate. Numerical implement shows the effectiveness of the proposed method in Section 4.


## 2. Regularization and Error Estimate

Letg^(ξ) denote the Fourier transform of the function g(y), which is defined as
(2)g^(ξ)=12π∫-∞∞e-iyξg(y)dy‍,i=-1.
The functions φ(y),φδ(y)∈L2(ℝ) are the exact and measured data for problem (1), respectively, and satisfy
(3)∥φδ(·)-φ(·)∥≤δ,
where ∥·∥ denotes the L2-norm and the constant δ>0 is the noise level. Assume that u(x,·)∈L2(ℝ) for all 0≤x<1 and there is the followinga priori bound;
(4)∥u(1,·)∥≤E,
where E is a positive constant.It is easy to know that for problem (1),
(5)u^(x,ξ)=cosh⁡(x|ξ|2-k2)φ^(ξ),
and equivalently,
(6)u(x,y)=12π∫-∞∞eiyξcosh⁡(x|ξ|2-k2)φ^(ξ)dξ‍.Note that the factorcosh⁡(x|ξ|2-k2) increases exponentially for 0<x<1 as |ξ|→+∞; a small distribution for the data φ(x) will be amplified infinitely by this factor and lead to the integral (6) blow-up. Therefore, recovering the temperature u(x,y) from the measured data φδ(x) is severely ill-posed.For simplicity [15], we decompose ℝ into the following parts I and W, where:
(7)I:={ξ∈ℝ,|ξ|≥k},W:={ξ∈ℝ,|ξ|≤k};
then L2(ℝ)=L2(I)⊕L2(W).Forξ∈W, we can take the regularization approximation solution in the frequency domain as
(8)u^mδ(x,ξ)=cos⁡(xk2-|ξ|2)φ^δ(ξ).Forξ∈I, we introduce an iteration scheme with the following form:
(9)u^mδ(x,ξ)=(1-λ)u^m-1δ(x,ξ)+λcosh⁡⁡(x|ξ|2-k2)φ^δ(ξ),hhhhhhhhhm=1,2,…,
where λ=λ(ξ)=e-|ξ|2-k2<1 plays an important role in the convergence proof; the initial guess is u^0δ(x,ξ). By using an elementary calculation for (9), we obtain
(10)u^mδ(x,ξ)=(1-λ)mu^0δ(x,ξ)+∑i=0m-1(1-λ)iλcosh⁡(x|ξ|2-k2)φ^δ(ξ)‍=(1-λ)mu^0δ(x,ξ)+(1-(1-λ)m)×cosh⁡(x|ξ|2-k2)φ^δ(ξ).
Therefore, the approximate solution of problem (1) has the following form in the frequency domain:
(11)u^mδ(x,ξ)={cos⁡(xk2-|ξ|2)φ^δ(ξ),|ξ|≤k,(1-λ)mu^0δ(x,ξ)+(1-(1-λ)m)×cosh⁡(x|ξ|2-k2)φ^δ(ξ),|ξ|≥k,
or equivalently,
(12)umδ(x,y)=12π∫-∞∞eixξu^mδ(x,ξ)dξ‍,
where u^mδ(x,ξ) is given by (11).Lemma 1 (see [16]).
For0<λ≤1 and m≥1, the following inequalities hold:
(13)(1-λ)mλ≤1m+1,1-(1-λ)mλ≤m.Lemma 2.
For0≤λ≤1,0≤α≤1 and m≥1, Lemma 1 can be strengthened as the following inequalities:
(14)(1-λ)mλα≤(m+1)-α,(15)1-(1-λ)mλα≤mα.Proof.
In fact, using the established results (13), we can get
(16)(1-λ)mλα≤[(1-λ)mλ]α≤(m+1)-α,1-(1-λ)mλα≤[1-(1-λ)mλ]α≤mα.Theorem 3.
Letu(x,y) be the exact solution of problem (1) and umδ(x,y) be its regularized approximation given by (12) with u0δ(x,y)=0. Assumptions (3) and (4) are satisfied and one chooses m=⌊E/δ⌋, where ⌊t⌋ denotes the largest integer not exceeding t; then there holds the following estimate:
(17)∥u(x,·)-umδ(x,·)∥≤3Exδ1-x+δ,0<x<1.Proof.
Due to the Parseval formula and the triangle inequality, we have(18)∥u(x,·)-umδ(x,·)∥2=∥u^(x,·)-u^mδ(x,·)∥2=∥u^(x,·)-u^mδ(x,·)∥L2(W)2+∥u^(x,·)-u^mδ(x,·)∥L2(I)2.
Case  1. While ξ∈W, combining (3), (6), and (12), we have
(19)∥u^(x,·)-u^mδ(x,·)∥L2(W)=∥cos⁡(xk2-|ξ|2)φ^(ξ)hhlh-cos⁡(xk2-|ξ|2)φ^δ(ξ)∥L2(W)≤∥φ^(ξ)-φ^δ(ξ)∥L2(W)≤∥φ^(ξ)-φ^δ(ξ)∥≤δ.
Case  2. For ξ∈I, combining (4), (6), (12), (14), and (15), we have
(20)∥u^(x,·)-u^mδ(x,·)∥L2(I)=∥cosh⁡(x|ξ|2-k2)φ^hhhhh-cosh⁡⁡(x|ξ|2-k2)(1-(1-λ)m)φ^δ∥≤∥cosh⁡(x|ξ|2-k2)(1-λ)mφ^∥+∥cosh⁡(x|ξ|2-k2)(1-(1-λ)m)(φ^-φ^δ)∥=∥cosh⁡(x|ξ|2-k2)(1-λ)mcosh⁡⁡(|ξ|2-k2)cosh⁡(|ξ|2-k2)φ^∥+∥cosh⁡(x|ξ|2-k2)(1-(1-λ)m)(φ^δ-φ^)∥≤Esup⁡ξ∈ℝcosh⁡(x|ξ|2-k2)(1-λ)mcosh⁡(|ξ|2-k2)+δsup⁡ξ∈ℝcosh⁡(x|ξ|2-k2)(1-(1-λ)m)≤2Esup⁡ξ∈ℝ(1-λ)me(x-1)|ξ|2-k2+δsup⁡ξ∈ℝ(1-(1-λ)m)ex|ξ|2-k2=2Esup⁡0<λ<1λ(1-x)(1-λ)m+δsup⁡0<λ<1λ-x(1-(1-λ)m)≤2E(m+1)x-1+δmx.
Due tom=⌊E/δ⌋, then m≤E/δ and m+1≥E/δ, therefore,
(21)∥u^(x,·)-u^mδ(x,·)∥L2(I)≤2E(Eδ)x-1+δ(Eδ)x=3Exδ1-x.
Combing inequalities (18), (19), and (21), the proof of this theorem is completed.Remark 4.
Obviously, Theorem3 could only solve the problem with the case 0<x<1. The stronger smoothness assumption of ∥u(1,·)∥ may obtain convergence rates for the endpoint x=1; see, for example, [10–12], and we omit the further discussions.


## 3. The Discrepancy Principle

In this section, we discuss ana posteriori stopping rule for iterative scheme (9) which is based on the discrepancy principle of Morozov [17, 18] in the following form:
(22)∥φδ-umδ(0,·)∥≤τδ≤∥φδ-um-1δ(0,·)∥,
where τ>1 is a constant and m denotes the regularization parameter. In the numerical experiments, we can take the iteration depth m which satisfies (22) firstly.Ifu^0δ(x,ξ)=0, thus (22) can be simplified to
(23)∥(1-λ)mφ^δ∥L2(I)≤τδ≤∥(1-λ)m-1φ^δ∥L2(I).Lemma 5.
The following inequality holds:(24)m≤2E(τ-1)δ.Proof.
Due to (4) and (14), we know
(25)τδ≤∥φδ-um-1δ(0,·)∥=∥(1-λ)m-1φ^δ∥L2(I)≤∥(1-λ)m-1(φ^δ-φ^)∥L2(I)+∥(1-λ)m-1φ^∥L2(I)≤∥φ^δ-φ^∥L2(I)+∥(1-λ)m-1cosh⁡(|ξ|2-k2)φ^cosh⁡(|ξ|2-k2)∥L2(I)≤δ+2Esup⁡0<λ<1(1-λ)m-1λ≤δ+2Em,
therefore,
(26)m≤2E(τ-1)δ.Lemma 6.
Settingωm(x,·)=u(x,·)-um(x,·), then the following inequality holds:
(27)∥ωm(x,·)∥L2(I)≤2∥ωm(1,·)∥L2(I)x∥ωm(0,·)∥L2(I)1-x.Proof.
Definingum(x,·)=(1-(1-λ)m)cosh⁡(x|ξ|2-k2)φ^(ξ), then we have
(28)ωm(0,·)=(1-λ)mφ^,ωm(1,·)=(1-λ)mcosh⁡(|ξ|2-k2)φ^.∥ωm(x,·)∥L2(I)2=∥(1-λ)mcosh⁡(x|ξ|2-k2)φ^∥L2(I)2=∫ξ∈I(1-e-|ξ|2-k2)2m(cosh⁡(x|ξ|2-k2)φ)2dξ‍≤∫ξ∈I(1-e-|ξ|2-k2)2m(ex|ξ|2-k2φ)2dξ‍=∫ξ∈I‍(1-e|ξ|2-k2)2mx(e|ξ|2-k2φ)2xHH×(1-e|ξ|2-k2)2m(1-x)(e|ξ|2-k2φ)2(1-x)dξ≤(∫ξ∈I((1-e-|ξ|2-k2)me|ξ|2-k2φ)2dξ‍)x×(∫ξ∈I((1-e-|ξ|2-k2)mφ)2dξ‍)1-x≤4(∫ξ∈I((1-e-|ξ|2-k2)mcosh⁡(|ξ|2-k2)φ)2dξ‍)x×(∫ξ∈I((1-e-|ξ|2-k2)mφ)2dξ‍)1-x=4∥ωm(1,·)∥L2(I)2x∥ωm(0,·)∥L2(I)2(1-x).Lemma 7.
The following inequality holds:(29)∥um(x,·)-u(x,·)∥L2(I)≤2(τ+1)1-xExδ1-x.Proof.
Due to (3), (4), and (27), we know that
(30)∥um(1,·)-u(1,·)∥L2(I)=∥(1-λ)mcosh⁡(|ξ|2-k2)φ∥≤E,∥um(0,·)-u(0,·)∥L2(I)=∥(1-λ)mφ∥L2(I)≤∥(1-e-|ξ|2-k2)m(φ-φδ)∥L2(I)+∥(1-e-|ξ|2-k2)mφδ∥L2(I)≤(τ+1)δ.
Combining (30), we have
(31)∥um(x,·)-u(x,·)∥L2(I)≤2(τ+1)1-xExδ1-x.Theorem 8.
Letu(x,y) be the exact solution of problem (1) and umδ(x,y) be its regularization approximation defined by (12) with u0δ(x,y)=0. If thea prioribound (4) is valid and the iteration (9) is stopped by the discrepancy principle (22), then
(32)∥umδ(x,·)-u(x,·)∥≤CEy/y0δ1-(y/y0)+δ,
where C=(2/(τ-1))x+2(τ+1)1-x.Proof.
According to the triangle inequality, (19), (24), and (29), we obtain that
(33)∥umδ(x,·)-u(x,·)∥≤∥umδ(x,·)-u(x,·)∥L2(I)+∥umδ(x,·)-u(x,·)∥L2(W)≤∥umδ(x,·)-um(x,·)∥L2(I)+∥um(x,·)-u(x,·)∥L2(I)+δ≤δmx+2(τ+1)1-xExδ1-x+δ≤δ(2E(τ-1)δ)x+2(τ+1)1-xExδ1-x+δ=((2τ-1)x+2(τ+1)1-x)Exδ1-x+δ.


## 4. Numerical Test

In this section, a simple numerical example is devised to verify the validity of the proposed method. We use the discrete Fourier transform and inverse Fourier transform (or FFT and IFFT algorithms) to complete our numerical experiment. We fix the intervala≤y≤b,  N denotes the number of discrete points.For an exact data functionφ(y), its discrete noisy version is
(34)φϵ=φ+ϵrandn(size(φ)),
where
(35)φ=(φ(y1),…,φ(yN)),yj=a+(b-a)(j-1)N+1,j=1,2…,N,δ=∥φϵ-φ∥l2:=1N∑i=1N|φϵ(yi)-φ(yi)|2‍.
The function “randn(·)” generates arrays of random numbers whose elements are normally distributed with mean 0, variance σ2=1. The absolute error ea(u) and the relative error er(u) are defined by
(36)ea(u):=∥umδ(x,·)-u(x,·)∥l2,(37)er(u):=∥umδ(x,·)-u(x,·)∥l2∥u(x,·)∥l2,
respectively.In the numerical experiment, we compute the approximationumδ(x,y) according to Theorem 3. And we can take the discrete points N=100, the number of wave k=1,a priori bound E=∥u(1,·)∥l2(ℝ)≈0.7, anda priori parameter m=⌊E/δ⌋. Thea posteriori parameter m was chosen according to formula (22) and τ=2 for calculation. Meanwhile, we take a=-6 and b=6 in the first numerical example. For Example 2, we take a=-4 and b=4.Example 1.
If we take the functionφ(y)=e-y2∈𝒮(ℝ), where 𝒮(ℝ) denotes the Schwartz function space, φ^(ξ)∈𝒮(ℝ) decays rapidly and formula (6) can be used to calculate u(x,y) with exact data directly. To observe the effect on different noise levels ϵ, we only take the case of k=1 at x=0.9.Table1 shows the comparison of the errors between the exact and regularization solutions for different ϵ, from which we can see that the smaller the ϵ is, the better the computed approximation is.Table 1
The errors between the exact and approximate solutions of Example1, with k=1 at x=0.9.
ϵ
1
e
-
4
1
e
-
3
1
e
-
2
1
e
-
1
m
4909
534
51
4
e
a
(
u
)
0.4435
0.5194
0.8220
2.1846
e
r
(
u
)
0.1021
0.1196
0.1893
0.5031Figure1 is the comparison ofa priori anda posteriori parameter choice rules for the exact u(x,y) and the approximate solution ukδ(x,y) at x=0.1,0.3,0.6,0.9 for the noise level ϵ=10-2. Here we also take the reasonablea priori bound E=∥u(1,·)∥l2(ℝ), and we can see that thea posteriori rule also works effectively.Example1. The regularization solution witha priori anda posteriori parameter choice rules for the noise level ϵ=10-2. (a) y=0.1, (b) y=0.3, (c) y=0.6, and (d) y=0.9, respectively.
(a)
(b)
(c)
(d)Example 2.
The function(38)u(x,y)=1π2sin(πky)cosh⁡(kπ2-1x)
is the exact solution of problem (1) with the Cauchy data u(0,y)=φ(y)=(1/π2)sin(πky) and ux(0,y)=0.Figure2 is the comparison of the exact solution u(x,y) and the approximation umδ(x,y) at different points x=0.1,0.3,0.6,0.9 and noise levels ϵ=10-2,10-3 for thea priori parameter choice rule. Here we take thea priori bound E=∥u(1,·)∥l2(ℝ), and the proposed method works well for thea priori parameter choice rule.Example2. The exact and regularized solutions for the different noise levels. (a) x=0.1, (b) x=0.3, (c) x=0.6, and (d) x=0.9  ϵ=10-2,10-3, respectively.
(a)
(b)
(c)
(d)Figure3 is the comparison of the differenta priori bound E=0.07,0.7,7 for the different points x=0.1,0.3,0.6,0.9 at the noise level ϵ=10-2. From this figure we can see that working with a wrong constant E would lead to a bad regularized solution. Therefore, a reasonablea priori bound E is very important for thea priori parameter choice rule.Example2. The exact and regularized solutions at (a) x=0.1, (b) x=0.3, (c) x=0.6, and (d) x=0.9 for the same noise level ϵ=10-2 but differenta priori bounds E=0.07,0.7,7, respectively.
(a)
(b)
(c)
(d)Figure4 is the comparison ofa priori anda posteriori parameter choice rules for the exact u(x,y) and the approximate solution umδ(x,y) at x=0.1,0.3,0.6,0.9 for the noise level ϵ=10-2. Here we take the reasonablea priori bound E=∥u(1,·)∥l2(ℝ) as previously mentioned, and thea posteriori rule also works effectively as expected.Example2. The regularization solution witha priori anda posteriori parameter choice rules for the noise level ϵ=10-2. (a) x=0.1, (b) x=0.3, (c) x=0.6, and (d) x=0.9, respectively.
(a)
(b)
(c)
(d)From Figures1–4, we concluded that the smaller the ϵ is, the better the computed approximation is, and the bigger the x is, the worse the computed approximation is. In addition, thea priori bound E has great influence on the numerical results. Although thea posteriori regularization parameter selection rule does not rely ona priori bound E, it also works well.


## 5. Conclusion

In this paper an iteration regularization method is given for solving the numerical analytic continuation problem on a strip domain. Thea priori anda posteriori rules for choosing a regularization parameter with strict theory analysis are presented. In numerical aspect, the comparison with different parameter choice rules shows that the proposed method works effectively.



---

*Source: 917972-2014-03-25.xml*