- Email: [email protected]
The Fourier–Bessel method for solving the Cauchy problem connected with the Helmholtz equation
Journal of Computational and Applied Mathematics 311 (2017) 183–193
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
The Fourier–Bessel method for solving the Cauchy problem connected with the Helmholtz equation Minghui Liu a,∗ , Deyue Zhang b , Xu Zhou c , Feng Liu a a
People’s Public Security University of China, China
b
School of Mathematics, Jilin University, China
c
College of Computer Science and Technology, Jilin University, China
article
abstract
info
Article history: Received 11 August 2014 Received in revised form 8 April 2016 Keywords: Cauchy problem Helmholtz equation Fourier–Bessel expansion Truncated SVD No meshing
This paper is concerned with the Cauchy problem connected with the Helmholtz equation. We investigate the denseness of Fourier–Bessel functions and propose a numerical method for approximating the solution to the problem. The convergence and stability are analyzed with truncated singular value decomposition method. Numerical experiments are also presented to show the effectiveness of our method. © 2016 Elsevier B.V. All rights reserved.
1. Introduction This paper is concerned with the Cauchy problem for the Helmholtz equation, which arises in many areas of science and engineering [1–5], such as wave propagation, vibration, electromagnetic scattering and so on. The direct problems, such as the Dirichlet, Neumann or the mixed boundary value problems connected with the Helmholtz equation, have been studied thoroughly in the past years. However, in some practical problems, only some noisy data on a part of the boundary instead of that on the whole boundary can be obtained, which leads to some inverse problems. It is well known that the Cauchy problem for the Helmholtz equation is an inverse problem and severely ill-posed [6–8], i.e. the solutions do not depend continuously on Cauchy data and a small perturbation in the data may cause large deviation in the solution. There are some numerical methods in the literature for solving the Cauchy problem connected with the Helmholtz equation. We refer the reader to [4,9–14] for the alternating iterative boundary element method, the conjugate gradient boundary element method, the boundary knot method, the fundamental solution method, the moment method, and homotopy perturbation method, etc. The Herglotz wavefunction method can be found in [15]. The main purpose of this paper is to provide a simple and effective numerical method for solving Cauchy problems connected with the Helmholtz equation. The main idea is to approximate the solution of the Cauchy problem by the Fourier–Bessel functions, which are solutions of the Helmholtz equation
1uN + k2 uN = 0,
in R2
(1.1)
of the form uN =
N
cl Jl (kr )eilθ ,
l=−N
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (M. Liu).
http://dx.doi.org/10.1016/j.cam.2016.07.023 0377-0427/© 2016 Elsevier B.V. All rights reserved.
(1.2)
184
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
where Jl (r ) is Bessel function of order l. Then we derive two equations for the coefficients on the specified boundary, which can be solved by a regularization method. Then the data for the solution on the unspecified boundary can be obtained easily. Roughly speaking, the main advantage of the Fourier–Bessel method over standard Galerkin methods is that they allow one to choose element sizes which are many wavelengths in diameter, while standard Galerkin methods with linear elements need of order 10 elements per wavelength to obtain reasonable accuracy. Thus at high frequency the total number of degrees of freedom is much reduced, and due to their high-order convergence, the high accuracy is obtained. The paper is organized as follows. In Section 2, we present the Fourier–Bessel function expansion approximation of the solution to the Cauchy problem, and obtain two equations for the coefficients. In Section 3, we solve the equations by the truncated singular value decomposition (T-SVD), and analyze the convergence and stability. Finally, several numerical examples are presented to show the effectiveness of our method. 2. Fourier–Bessel function approximation 2.1. Formulation of the problem Let D ⊂ R2 be a simply connected and bounded domain with a regular boundary ∂ D, and Γ be a portion of the boundary ∂ D. Consider the following Cauchy problem: given Cauchy data fD and fN on Γ , find u on Σ = ∂ D \ Γ such that u satisfies the Helmholtz equation:
1u + k2 u = 0, in D, u = fD , on Γ ,
∂u = fN , ∂n
(2.1) (2.2)
on Γ ,
(2.3)
where n is the unit normal with respect to the boundary ∂ D, and k > 0 is the wave number. In this paper, we assume that −k2 is not an eigenvalue of the Laplacian operator with the homogeneous Neumann boundary condition. Without loss of generality, we assume that the measured data fD ∈ H 1 (Γ ) and fN ∈ L2 (Γ ), and the Cauchy problem has a unique solution u in H 3/2 (D). The approximation can be split into two steps. First, we approximate the solution of the Cauchy problem by a single layer potential. Second, we approximate the single layer potential by Bessel functions. 2.2. Single layer potential approximation Define the single layer potential vϕ by
vϕ (x) :=
∂B
Φ (x, y )ϕ(y )ds(y ),
x ∈ B,
where B is a Lipschitz domain, vϕ ∈ L2 (∂ B) and Φ (x, y ) = equation. Then we have (see [16]).
(2.4) i (1) H 4 0
(k|x − y |) is the fundamental solution of the Helmholtz
Lemma 2.1. Let u ∈ H 3/2 (D) satisfy the Helmholtz equation. Let B be a bounded and simply connected domain with ∂ B ∈ C 2 , ¯ ⊂⊂ B. Then for every ε > 0, there exists a single layer potential of the form and D
vϕ (x) :=
∂B
Φ (x, y )ϕ(y )ds(y ),
x ∈ B,
(2.5)
for some ϕ ∈ L2 (∂ B), such that
∥vϕ − u∥H 3/2 (D) 6 ε,
(2.6)
∂vϕ ∂u ∥vϕ − u∥L2 (∂ D) + − 6 ε. ∂n ∂ n L2 (∂ D)
(2.7)
and
∂v
Proof. It is sufficient to prove that the set {( ∂ nϕ + ivϕ )|∂ D : ϕ ∈ L2 (∂ B)} is dense in L2 (∂ D). Let φ ∈ L2 (∂ D), satisfying
∂D
∂vϕ + ivϕ φ¯ ds(x) = 0, ∂n
for all ϕ ∈ L2 (∂ B).
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
185
Then
∂ Φ (x, y )ϕ(y )ds(y ) + i Φ (x, y )ϕ(y )ds(y ) φ(x)ds(x) ∂ D ∂ nx ∂ B ∂B ∂ Φ (x, y ) + iΦ (x, y ) φ(x)ds(x)ds(y ) = 0, = ϕ(y ) ∂ nx ∂B ∂D
for all ϕ ∈ L2 (∂ B). Define
w( ˜ y ) :=
∂D
∂ Φ (x, y ) + iΦ (x, y ) φ(x)ds(x), ∂ nx
y ∈ ∂ B,
then it is readily seen that w( ˜ y ) = 0 on ∂ B. This means that w ˜ solves the boundary value problem
1w ˜ + k2 w ˜ = 0 in R2 \ D¯ ⊃ R2 \ B¯ , w ˜ = 0 on ∂ B.
¯ And a unique continuation argument Since w ˜ satisfies Sommerfeld radiation condition, it is obvious that w( ˜ y ) = 0 in R2 \ B. 2 ¯ yields that w( ˜ y ) = 0 in R \ D. From the jump relation, we have ¯ w ˜+ −w ˜ − = φ,
∂w ˜+ ∂w ˜+ ¯ − = −iφ, ∂n ∂n
∂w ˜
on ∂ D.
∂w ˜
Since ∂ n+ + iw ˜ + = 0 on ∂ D, we have that ∂ n− + iw ˜ − = 0 on ∂ D. Therefore, w ˜ = 0 in D, and φ = 0 on ∂ D. From the denseness, we have that for every ε > 0, there exists a single layer potential vϕ of the form (2.4) for some ϕ ∈ L2 (∂ B), such that
∂vϕ ∂u 6 ε. + i v − + iu ϕ ∂n 2 ∂n L (∂ D) From the interior regularity results of the Helmholtz equation, we have
∥vϕ − u∥H 3/2 (D) 6 ε. By using the trace theorem, we obtain the estimate (2.7).
2.3. Bessel function approximation Now we approximate vϕ in D by Bessel functions. First, we introduce the Vekua operator [17–19] in the case of Helmholtz equation 1u + k2 u = 0. Let Ω ∈ C, Ω ∗ = {z ∈ C : z¯ ∈ Ω },
G(t , t ∗ , z , z ∗ ) := J0 (k (z − t )(z ∗ − t ∗ )),
∂ G(t , t ∗ , z , z ∗ ). ∂t Let v be an analytic function in Ω , then for fixed z0 ∈ Ω , we have the Vekua operator defined as z ReV[v, z0 ](z ) := Re{G(z , z¯0 , z , z¯ )v(z ) + v(t )H (z , z¯0 , z , z¯ )dt }, z = x1 + ix2 ∈ Ω . H (t , t ∗ , z , z ∗ ) := −
(2.8)
z0
And we have the following results (see [18]). Lemma 2.2 (Vekua). Let Ω be simply connected and fixed z0 ∈ Ω . Then there exists a unique analytic function v in Ω with v(z0 ) real, such that u(x1 , x2 ) = ReV[v, z0 ](z ),
z = x1 + ix2 ∈ Ω .
If Ω is bounded, it is readily seen that
∥ReV[v, z0 ](z )∥L∞ (Ω ) 6 ∥G∥∞ ∥v∥L∞ (Ω ) +
z
∥v∥L∞ (Ω ) ∥H ∥∞ d|t | 6 KV ∥v∥L∞ (Ω ) , z0
where ∥G∥∞ , ∥H ∥∞ are the supreme of these functions on Ω × Ω ∗ × Ω × Ω ∗ , and KV only depends on Ω and k.
186
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
Besides, a direct calculation gives that ReV[z n , 0](z ) =
n!2n
cos(nθ )Jn (kr ), kn n!2n ReV[iz n , 0](z ) = − n sin(nθ )Jn (kr ), k with x1 = r cos(θ ), x2 = r sin(θ ). Now, we give the approximation theorem in L∞ with Fourier–Bessel expansion. Theorem 2.1. Assume that u is the solution of the Helmholtz equation (2.1) in Bρ := {x ∈ R2 : |x| < ρ}, then there exists a Fourier–Bessel expansion of the form uN (r , θ ) =
N
cl Jl (kr )eilθ ,
(2.9)
l=−N
and we have the estimate
∥u − uN ∥L∞ (Br ) 6 C1 τ −N ,
(2.10)
∥∇(u − uN )∥L∞ (Br ) 6 C2 τ −N ,
(2.11)
where r < R < ρ , τ = R/r > 1, and C1 , C2 are constants. Proof. By Lemma 2.2, there exists a unique analytic function v ∈ Bρ , such that u = ReV[v, 0](z ). Then there exists polynomials pN of degree N, such that [15,19,18]
∥v − pN ∥L∞ (Br ) 6 c˜1 r −N ∥v∥L∞ (∂ Bρ ) ,
(2.12)
∥∇(v − pN )∥L∞ (Br ) 6 c˜2 r −N ∥v∥L∞ (∂ Bρ ) ,
(2.13)
where c˜1 and c˜2 only depends on r and ρ . Then by using the equality J−l (kr ) = (−1)l Jl (kr ), l = 1, 2, . . . , we have
∥u − uN ∥L∞ (Br )
N ilθ cl Jl (kr )e = u − ∞ l=−N L (Br ) N N (2) (1) al Jl (kr ) sin(lθ ) al Jl (kr ) cos(lθ ) + = u − ∞ l=1 l =0 L (Br ) N N (2) (1) a˜ l iz l a˜ l z l + 6 KV v − ∞ l =1 l=0 L
6 c˜1 KV τ (1)
(1)
−N
cl +(−1)l c−l
∥v∥L∞ (∂ Bρ ) 6 C1 τ (2)
cl −(−1)l c−l
, al = where a0 = c0 , al = 2 2i Similarly, we can obtain the estimate (2.11).
−N
,
(1)
, a˜ l
(Br )
=
kl (1) a , 2l l! l
(2)
a˜ l
= − 2kl l! a(l 2) , l = 1, . . . , N. l
Thus, we obtain the main result in this section. Theorem 2.2. Let u ∈ H 3/2 (D) satisfies Helmholtz equation. Then for every ε > 0, there exists a Fourier–Bessel expansion uN of the form (2.9), such that
∥u − uN ∥H 3/2 (D) 6 ε,
(2.14)
∂ u ∂ uN − ∥u − uN ∥L2 (∂ D) + 6 ε. ∂n ∂ n L2 (∂ D)
(2.15)
and
¯ ⊂⊂ Br ⊂ Bρ ⊂ B, where r < R < ρ . Proof. Choose B a bounded and simply connected domain with ∂ B ∈ C 2 , such that D From Lemma 2.1, we have that for every ε > 0, there exists a single layer potential vϕ of the form (2.4) for some ϕ ∈ L2 (∂ B), such that ∥vϕ − u∥H 3/2 (D) 6 ε.
(2.16)
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
187
By using Theorem 2.1 for vϕ in Bρ , then for ε > 0, there exists a Fourier–Bessel expansion uN of the form (2.9), such that
∥vϕ − uN ∥L∞ (Br ) 6 C1 τ −N 6 ε, ∥∇(vϕ − uN )∥L∞ (Br ) 6 C2 τ −N 6 ε, where τ = R/r > 1. Then it is readily seen that
∥vϕ − uN ∥L2 (∂ D) 6 C˜ 1 ε, and
∂vϕ − uN 6 C˜ 2 ε. ∂n 2 L (∂ D) From the interior regularity of the Helmholtz equation, we have that
∥vϕ − uN ∥H 3/2 (D) 6 C3 ε.
(2.17)
Then by (2.16) and (2.17),
∥u − uN ∥H 3/2 (D) 6 ∥u − vϕ ∥H 3/2 (D) + ∥vϕ − uN ∥H 3/2 (D) 6 C ε. This means (2.14) is valid. By using the trace theorem, the estimate (2.15) is obtained.
Now we derive the numerical approximation of uN . First, define the trace operator An : C2N +1 → L2 (Γ ) × L2 (Γ ) by
N
l=−N
An cn :=
cl Jl (kr )eilθ
N ∂ cl Jl (kr )e ∂ n l=−N
, ilθ
x ∈ Γ,
(2.18)
where cn = (c−N , c−N +1 , . . . , cN −1 , cN ) ∈ C2N +1 . Then we have the following theorem. Theorem 2.3. The operator An : C2N +1 → L2 (Γ ) × L2 (Γ ) defined by (2.18) is injective and compact. Proof. Let An cn = 0, then there exists a Fourier–Bessel expansion uN (x) =
N
∂ uN | ∂n Γ
l=−N
cl Jl (kr )eilθ , such that uN |Γ = 0 and
= 0. From the Green formula and the analyticity of the harmonic polynomial, it is readily seen that uN = 0 in R2 , and so cn = 0. Therefore, the operator An is injective. 2N +1 Let {c (k) }∞ , then there exists a convergent subsequence {c (kj ) }∞ k=1 be a bounded sequence in C j=1 , such that limj→∞ |c (kj ) − c˜ | = 0 for some c˜ ∈ C2N +1 .
√
Using the well known results of Bessel functions: |J0 (kr )| < 1, |Jl (kr )| < 1/ 2 for l = 1, 2, . . . , and J−l (kr ) = (−1)l Jl (kr ) for l = 1, 2, . . . , we have
∥J0 (kr )eilθ ∥2L2 (∂ D) 6 2π ,
∥Jl (kr )eilθ ∥2L2 (∂ D) 6 π ,
l = ±1, ±2, . . . ,
and realizing that Jl′ (r ) = Jl−1 (r ) + Jl+1 (r ), we have
2 ∂ Jl (kr )eilθ ∂n 2
L (∂ D)
2 ∂ ilθ 6 C Jl (kr )e ∂r 2
L (∂ D)
2 ∂ ilθ + Jl (kr )e ∂θ 2
L (∂ D)
6C
k2 2
+
l2
2
,
where C is a constant only depending on ∂ D. So by using the Cauchy inequality, we have
∥An c
(kj )
− An c˜ ∥
2 L2 (Γ )×L2 (Γ )
2 2 N N ∂ (kj ) (kj ) ilθ ilθ = (cl − c˜l )Jl (kr )e + (cl − c˜l )Jl (kr )e l=−N 2 ∂ n l=−N 2 L (Γ ) L (Γ ) 2 N ∂ ilθ 6 |c (kj ) − c˜ |2 ∥Jl (kr )eilθ ∥2L2 (∂ D) + ∂ n Jl (kr )e 2 L (∂ D) l=−N 2 k N ( N + 1)(2N + 1) (1 + 2N ) + . 6 |c (kj ) − c˜ |2 2π (1 + N ) + C 2
So limj→∞ ∥An c
(kj )
− An c˜ ∥L2 (Γ )×L2 (Γ ) = 0. Therefore, An is compact.
6
188
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
3. A numerical method for solving the problem Numerically, we should solve the coefficients by the equation
(An cn )(x) = b(x),
x ∈ Γ,
(3.1)
where b(x) = (f (x), g (x))T ∈ L2 (Γ ) × L2 (Γ ). But (3.1) is a severely ill-posed problem and we need a regularization method. In practice, the measured Cauchy data b(x) is often perturbed, that is, the measured data is bδ with ∥bδ − b∥ 6 δ . In this section, we introduce the truncated singular value decomposition (T-SVD) method to overcome the ill-posedness of (3.1) with perturbed measured data, i.e., we solve the equation An cn = bδ .
(3.2)
First, we introduce the singular system of An : C → L (Γ ) × L (Γ ). Let µ1 > µ2 > · · · µ2N +1 > 0 be the singular +1 +1 values of An , and (xj )2N ⊂ C2N +1 , (yj )2N ⊂ L2 (Γ ) × L2 (Γ ) be orthonormal systems with j =1 j =1 2N +1
An xj = µj yj ,
A∗n yj = µj xj ,
2
2
j = 1, . . . , 2N + 1,
where A∗n : L2 (Γ ) × L2 (Γ ) → C2N +1 is the adjoint operator of An . Then for a chosen α = α(δ) > 0, cnα(δ),δ =
1 (bδ , yj )xj µj 2
(3.3)
µj >α
is the T-SVD solution of (3.2). And we have the following theorem. Theorem 3.1. (a) There exists a function h ∈ L2 (Γ ) × L2 (Γ ), such that for every ε > 0,
|cn − A∗n h| < ε. (b) Let ∥h∥L2 (Γ )×L2 (Γ ) 6 E, then for the choice α(δ) = c δ/E for some c > 0, we have √ √ 1 δ E + ε. |cnα(δ),δ − cn | 6 √ + c
(3.4)
c
Proof. (a) Since An : C2N +1 → L2 (Γ ) × L2 (Γ ) is injective and compact, the range of A∗n is dense in C2N +1 . Therefore, for every ε > 0, there exists a function h ∈ L2 (Γ ) × L2 (Γ ), such that
|cn − A∗n h| < ε. (b) The proof can be obtained by the combination of (a) and Theorem 2.9 in book [20]. δ
δ
Numerically, the measured data b is in the discretization form (f (x1 ), f (x2 ), . . . , f δ (xm ), g δ (x1 ), g δ (x2 ), . . . , g δ (xm ))T := β δ , with Γm := {x1 , x2 , . . . , xm } ⊂ Γ , and m large enough to approximate the space L2 (Γ ) × L2 (Γ ). Denote Qm : L2 (Γ ) × L2 (Γ ) → Cm mapping bδ to β δ , then the problem (3.2) is solved in the finite space in the form Qm An cnδ = β δ ,
δ
(3.5)
δ
with |β − β| 6 δ , where β is the solution of (3.5) without noise. Now we apply the T-SVD method [21] to the discretization case (3.5). Note that Qm An : C2N +1 → Cm is a matrix, and denote it by A. We will show that a classical T-SVD method for matrix will do for (3.5) and that the solution of (3.5) will converge to that of (3.2) when m → ∞. The well known SVD of A can be written as A = U
Σ 0
V ∗ with Um×(2N +1) and V(2N +1)×(2N +1) unitary. Denote U =
(u1 , u2 , . . . , um ), V = (v1 , v2 , . . . , v2N +1 ), Σ = diag(λ1 , λ2 , . . . , λ2N +1 ), and assume that λ1 > λ2 > · · · > λ2N +1 > 0. Let 1 −1 −1 ˜ = diag(λ− α(δ) := Λ2 > λ2K for chosen K (1 6 K 6 2N + 1), and Σ 1 , λ2 , . . . , λK , 0, . . . , 0). Then the T-SVD solution of (3.5) can be written as
˜ U ∗ β δ := A˜ β δ . cnα(δ),δ = VΣ ,m
(3.6)
And a similar error estimate to Theorem 3.1 can be obtained. Theorem 3.2. Let cn,m be the solution of (3.5) with β the righthand side, and there exists z ∈ Cm with |z | 6 E, such that cn,m = A∗ z, then
δ + ΛE . |cnα(δ),δ − cn,m | 6 ,m Λ √ If we choose Λ = c δ/E for some c > 0, then √ √ 1 α(δ),δ |cn,m − cn,m | 6 √ + c δE . c
(3.7)
(3.8)
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
189
Proof.
|cnα(δ),δ − cn,m | = |A˜ β δ − A−1 β| ,m = |A˜ β δ − A˜ β + A˜ β − A−1 β| ˜ β δ − A˜ β| + |A˜ β − A−1 β|. 6 |A It is obvious that |A˜ β δ − A˜ β| 6 Λ−1 δ , and we have that
+1 +1 2N 2N 1 1 −1 ˜ (β, uj )vj = (An cn,m , uj )vj |Aβ − A β| = j=K +1 λj j=K +1 λj +1 +1 2N 2N = (cn,m , vj )vj = (Az , vj )vj j=K +1 j=K +1 +1 +1 2N 2N = (z , A∗ vj )vj = (z , λj uj )vj j=K +1 j=K +1 +1 2N λj (z , uj )vj 6 λK +1 |z | 6 ΛE . = j=K +1 So,
δ + ΛE . Λ √ And if we choose Λ = c δ/E for some c > 0, it is readily seen that √ √ 1 + c δE . |cnα(δ),δ − c | 6 √ n ,m ,m |cnα(δ),δ − cn,m | 6 ,m
c
Remark 1. It is easy to implement the numerical method. It is nothing but the traditional SVD for the matrix A, which is the discretization form of the trace operator An . Finally, we have the convergence result. Lemma 3.1. Let the assumptions of Theorems 3.1 and 3.2 hold, then the solution of (3.5) converges to that of (3.2) with the same righthand side bδ and the same regularization parameter α(δ) when m tends to infinity, i.e. lim |cnα(δ),δ − cnα(δ),δ | = 0. ,m
m→∞
Proof. It is the direct result of Remark 2.3 in [22].
Then it is readily seen that Theorem 3.3. Let the assumptions of Theorems 3.1 and 3.2 hold, then for every ε > 0, choose Λ = m sufficiently large, we have
|cnα(δ),δ − cn | 6 ,m
1
√
√ + c
c
√ δ E + ε.
√
c δ/E for some c > 0, and
(3.9)
Since the approximation of the coefficients is obtained, the approximation of the solution of the Helmholtz equation can be written as uα(δ),δ n,m (x) =
N
α(δ),δ
cl,m
Jl (kr )eilθ ,
x ∈ R2 .
l=−N
And we have the main result in the paper.
(3.10)
190
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
(a) Re(uN ) and u.
(b) Re(∂n uN ) and ∂n u. Fig. 1. The exact solution u and the numerical solution uN on Σ with k = 5 for Example 1.
Theorem 3.4. Let the assumptions of Theorems 3.1 and 3.2 hold, then for every ε > 0, choose Λ = m sufficiently large, the following estimate on boundary Σ holds,
√
c δ/E for some c > 0, and
2 ∂ uα(δ),δ ∂u n ,m 6 (1/c + c + 2)δ E − − u∥ + ∂n ∂n 2 L (Σ ) 2 k N (N + 1)(2N + 1) · 2π (1 + N ) + C (1 + 2N ) + + ε,
∥uα(δ),δ n,m
2 L2 (Σ )
2
(3.11)
6
where C only depends on ∂ D. Proof. By using the Cauchy inequality, for any x ∈ ∂ D, we have
2 ∂ uα(δ),δ ∂u n ,m − u∥ − + ∂n ∂n 2 L (∂ D) N ∂(Jl (kr )eilθ ) 2 ilθ 2 α(δ),δ 2 ∥Jl (kr )e ∥L2 (∂ D) + 6 |cn,m − cn | 2 ∂n L (∂ D) l=−N 2 k N (N + 1)(2N + 1) 6 ((1/c + c + 2)δ E + ε) 2π (1 + N ) + C (1 + 2N ) + .
∥uα(δ),δ n,m
2 L2 (∂ D)
2
So, the estimate (3.11) is obvious.
6
4. Numerical examples In this section, we present two examples to demonstrate the competitiveness of our algorithm. Example 1. Consider the case in which the exact solution to Cauchy problem is u(x) =
1 2k2 2
√
sin( 2kx2 )(ekx1 + e−kx1 )
(see [12,15]). Let D = {(x1 , x2 )|(x1 − 0.5)2 + x22 < 0.52 , x2 > 0}, Γ = {(x1 , x2 )|(x1 − 0.5) + x22 = 0.52 , x2 > 0} and Σ = ∂ D \ Γ = {(x1 , x2 )|0 6 x1 6 1, x2 = 0}. Fig. 1 shows the exact solution and the numerical solution for various levels of noise with k = 5. Table 1 shows the relative L2 error of uN and ∂ uN on Σ , respectively. From the figure and the table, we can conclude that the numerical solution is a stable approximation of the exact solution, and it converges to the exact solution when the noise decreases. Fig. 2 shows the numerical solution for various levels of wavenumber with noise level σ = 0 and σ = 5%, respectively. The result tells us that the solution for smaller k can be more stable. Example 2. Consider the unit disc D = {(x1 , x2 )|x21 + x22 < 1}. Let Γ = {x ∈ ∂ D|0 < θ (x) < Θ } and Σ = ∂ D \ Γ = {x ∈ ∂ D|Θ < θ (x) < 2π }, where θ (x) is the polar angle of x and Θ is a specified angle. In this example, we observe the effect of Θ on the numerical solution. Choose u(x) = exp(ix · d ) a plane wave, where d = (cos(π /3), sin(π /3)).
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
191
Table 1 Errors for boundary data with different noise levels for Example 1. Noise level
0 1% 3% 5%
∥uN − u∥L2 (Σ ) k=1
k=3
k=5
∥∂ uN − ∂ u∥L2 (Σ ) /∥∂ u∥L2 (Σ ) k=1 k=3
k=5
3.11e−08 3.23e−03 9.19e−03 1.37e−02
1.71e−07 3.54e−03 1.20e−02 1.40e−02
2.70e−06 7.11e−03 2.44e−02 3.37e−02
1.91e−07 9.18e−03 2.19e−02 2.85e−02
6.06e−06 4.65e−03 1.77e−02 2.73e−02
(a) Re(uN ) and u, for σ = 0.
1.09e−06 1.48e−02 4.02e−02 5.65e−02
(b) Re(uN ) and u, for σ = 5%.
Fig. 2. The exact solution u and the numerical solution uN on Σ with σ = 0 and σ = 5% for Example 1.
(a) Re(uN ) and Re(u), for σ = 5% and k = 1.
(b) Re(∂n uN ) and Re(∂n u), for σ = 5% and k = 1.
Fig. 3. The exact solution u and the numerical solution uN on Σ with k = 1, 5% noise and different values of Θ for Example 2.
Figs. 3 and 4 compare the numerical results on the whole boundary ∂ D for different values of Θ , where the wavenumber k = 1 and the noise level σ = 5% in Fig. 3, and wavenumber k = 5 and the noise level σ = 1% in Fig. 4. From these figures, it is readily seen that large Θ makes the numerical reconstruction more stable and accurate. Tables 2 and 3 show the relative L2 error of uN and ∂ uN on Σ , respectively. The results tell us that larger Θ makes the reconstruction better. 5. Conclusion In this paper, we study the Fourier–Bessel function expansion method for solving the Cauchy problem connected with the Helmholtz equation. The numerical method is based on the truncated singular value decomposition (T-SVD). Convergence
192
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
(a) Re(uN ) and Re(u), for σ = 1% and k = 5.
(b) Re(∂n uN ) and Re(∂n u), for σ = 1% and k = 5.
Fig. 4. The exact solution u and the numerical solution uN on Σ with k = 5, 1% noise and different values of Θ for Example 2. Table 2 Relative L2 errors on boundary Σ for different values of Θ and noise levels for Example 2: ∥uN − u∥L2 (Σ ) /∥u∥L2 (Σ ) . Noise level
0 1% 3% 5%
k = 0.5
k=5
Θ = π/2
Θ=π
Θ = 3π/2
Θ = π/2
Θ=π
Θ = 3π/2
7.01e−08 2.81e−02 2.07e−02 2.46e−02
1.11e−09 2.03e−03 1.62e−02 8.64e−03
9.29e−12 5.45e−04 5.12e−03 3.35e−03
4.99e−03 2.75e−01 3.46e−01 4.22e−01
2.67e−06 8.62e−02 1.23e−01 1.29e−01
1.89e−09 1.54e−03 6.97e−03 1.27e−02
Table 3 Relative L2 errors on boundary Σ for different values of Θ and noise levels for Example 2: ∥∂ uN − ∂ u∥L2 (Σ ) /∥∂ u∥L2 (Σ ) . Noise level
0 1% 3% 5%
k = 0.5
k=5
Θ = π/2
Θ=π
Θ = 3π/2
Θ = π/2
Θ=π
Θ = 3π/2
8.68e−07 1.31e−01 9.23e−02 1.04e−01
1.40e−08 1.56e−02 7.38e−02 1.19e−02
1.96e−10 4.91e−03 5.20e−02 2.80e−02
8.80e−03 2.46e−01 3.28e−01 3.59e−01
6.60e−06 1.01e−01 1.17e−01 1.20e−01
8.62e−09 4.64e−03 1.61e−02 2.35e−02
and stability are analyzed. The effectiveness of the numerical solution depends on the domain, the wavenumber, and the number of Bessel functions used. The method does not need meshing and computes fast. Acknowledgments The first author was supported by the Fundamental Research Funds for the Central Universities Grant 2014JKF03104, and the second author was supported by the National Natural Science Foundation of China (NSFC) Grant 11271159. The authors appreciate very much for the reviewers’ valuable suggestions and comments. References [1] M.R. Bai, Application of bem (boundary element method)-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries, J. Acoust. Soc. Am. 92 (1992) 533–549. [2] W. Hall, X. Mao, A boundary element investigation of irregular frequencies in electromagnetic scattering, Eng. Anal. Bound. Elem. 16 (3) (1995) 245–252. [3] B.-K. Kim, J.-G. Ih, On the reconstruction of the vibro-acoustic field over the surface enclosing an interior space using the boundary element method, J. Acoust. Soc. Am. 100 (5) (1996) 3003–3016. [4] L. Marin, L. Elliott, P. Heggs, D. Ingham, D. Lesnic, X. Wen, An alternating iterative algorithm for the Cauchy problem associated to the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 192 (5) (2003) 709–722. [5] W. Chen, Z.-J. Fu, Boundary particle method for inverse Cauchy problems of inhomogeneous Helmholtz equations, J. Mar. Sci. Technol. 17 (3) (2009) 157–163. [6] E.R.G. Alessandrini, L. Rondi, S. Vessella, The stability for the Cauchy problem for elliptic equations, Inverse Problems 25 (3) (2009) 1–47. [7] V.I.T. Hrycak, Increased stability in the continuation of solutions to the Helmholtz equation, Inverse Problems 20 (3) (2004) 697–712.
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
193
[8] V. Isakov, Inverse Problems for Partial Differential Equations, in: Applied Mathematical Sciences, vol. 120, Springer-Verlag, New York, 1998, http://dx.doi.org/10.1007/978-1-4612-5338-9. [9] B. Jin, Y. Zheng, A meshless method for some inverse problems associated with the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 195 (19) (2006) 2270–2288. [10] L. Marin, D. Lesnic, The method of fundamental solutions for the Cauchy problem associated with two-dimensional Helmholtz-type equations, Comput. Struct. 83 (4) (2005) 267–278. [11] T. Wei, Y. Hon, L. Ling, Method of fundamental solutions with regularization techniques for Cauchy problems of elliptic operators, Eng. Anal. Bound. Elem. 31 (4) (2007) 373–385. [12] T. Wei, H. Qin, R. Shi, Numerical solution of an inverse 2D Cauchy problem connected with the Helmholtz equation, Inverse Problems 24 (3) (2008) 035003. [13] A. Zakeri, A. Aminataei, Q. Jannati, Application of He’s Homotopy perturbation method for Cauchy problem of ill-posed nonlinear diffusion equation, Discrete Dyn. Nat. Soc. (2010) http://dx.doi.org/10.1155/2010/780207. [14] M.A. Jafari, A. Aminataei, Exact solution for the inverse problem of finding a source term in a parabolic equation, J. Stat. Math. 2 (2) (2011) 37–39. [15] D. Zhang, F. Ma, E. Zheng, A Herglotz wave function method for the inverse Cauchy problem of the Helmholtz equation, J. Comput. Appl. Math. 237 (1) (2013) 215–222. [16] Y. Sun, D. Zhang, F. Ma, A potential function method for the Cauchy problem of elliptic operators, J. Math. Anal. Appl. 395 (1) (2012) 164–174. [17] I. Vekua, New Methods of Solution of Elliptic Equations, Gostekhizdat, Moscow, 1948. [18] T. Betcke, A gsvd formulation of a domain decomposition method forplanar eigenvalue problems, IMA J. Numer. Anal. 27 (3) (2007) 451–478. [19] J. Melenk, Operator adapted spectral element methods i: harmonic and generalized harmonic polynomials, Numer. Math. 84 (1) (1999) 35–69. [20] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, in: Applied Mathematical Sciences, vol. 120, Springer-Verlag, New York, 1996, http://dx.doi.org/10.1007/978-1-4612-5338-9. [21] M. Liu, F. Ma, Numerical method for inverse scattering in two-layered background in near-field optics, Inverse Probl. Sci. Eng. 23 (1) (2015) 130–154. http://dx.doi.org/10.1080/17415977.2014.890611. [22] C. Vogel, Optimal choice of a truncation level for the truncated svd solution of linear first kind integral equations when data are noisy, SIAM J. Numer. Anal. 23 (1) (1986) 109–117.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
The Fourier–Bessel method for solving the Cauchy problem connected with the Helmholtz equation Minghui Liu a,∗ , Deyue Zhang b , Xu Zhou c , Feng Liu a a
People’s Public Security University of China, China
b
School of Mathematics, Jilin University, China
c
College of Computer Science and Technology, Jilin University, China
article
abstract
info
Article history: Received 11 August 2014 Received in revised form 8 April 2016 Keywords: Cauchy problem Helmholtz equation Fourier–Bessel expansion Truncated SVD No meshing
This paper is concerned with the Cauchy problem connected with the Helmholtz equation. We investigate the denseness of Fourier–Bessel functions and propose a numerical method for approximating the solution to the problem. The convergence and stability are analyzed with truncated singular value decomposition method. Numerical experiments are also presented to show the effectiveness of our method. © 2016 Elsevier B.V. All rights reserved.
1. Introduction This paper is concerned with the Cauchy problem for the Helmholtz equation, which arises in many areas of science and engineering [1–5], such as wave propagation, vibration, electromagnetic scattering and so on. The direct problems, such as the Dirichlet, Neumann or the mixed boundary value problems connected with the Helmholtz equation, have been studied thoroughly in the past years. However, in some practical problems, only some noisy data on a part of the boundary instead of that on the whole boundary can be obtained, which leads to some inverse problems. It is well known that the Cauchy problem for the Helmholtz equation is an inverse problem and severely ill-posed [6–8], i.e. the solutions do not depend continuously on Cauchy data and a small perturbation in the data may cause large deviation in the solution. There are some numerical methods in the literature for solving the Cauchy problem connected with the Helmholtz equation. We refer the reader to [4,9–14] for the alternating iterative boundary element method, the conjugate gradient boundary element method, the boundary knot method, the fundamental solution method, the moment method, and homotopy perturbation method, etc. The Herglotz wavefunction method can be found in [15]. The main purpose of this paper is to provide a simple and effective numerical method for solving Cauchy problems connected with the Helmholtz equation. The main idea is to approximate the solution of the Cauchy problem by the Fourier–Bessel functions, which are solutions of the Helmholtz equation
1uN + k2 uN = 0,
in R2
(1.1)
of the form uN =
N
cl Jl (kr )eilθ ,
l=−N
∗
Corresponding author. E-mail addresses: [email protected], [email protected] (M. Liu).
http://dx.doi.org/10.1016/j.cam.2016.07.023 0377-0427/© 2016 Elsevier B.V. All rights reserved.
(1.2)
184
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
where Jl (r ) is Bessel function of order l. Then we derive two equations for the coefficients on the specified boundary, which can be solved by a regularization method. Then the data for the solution on the unspecified boundary can be obtained easily. Roughly speaking, the main advantage of the Fourier–Bessel method over standard Galerkin methods is that they allow one to choose element sizes which are many wavelengths in diameter, while standard Galerkin methods with linear elements need of order 10 elements per wavelength to obtain reasonable accuracy. Thus at high frequency the total number of degrees of freedom is much reduced, and due to their high-order convergence, the high accuracy is obtained. The paper is organized as follows. In Section 2, we present the Fourier–Bessel function expansion approximation of the solution to the Cauchy problem, and obtain two equations for the coefficients. In Section 3, we solve the equations by the truncated singular value decomposition (T-SVD), and analyze the convergence and stability. Finally, several numerical examples are presented to show the effectiveness of our method. 2. Fourier–Bessel function approximation 2.1. Formulation of the problem Let D ⊂ R2 be a simply connected and bounded domain with a regular boundary ∂ D, and Γ be a portion of the boundary ∂ D. Consider the following Cauchy problem: given Cauchy data fD and fN on Γ , find u on Σ = ∂ D \ Γ such that u satisfies the Helmholtz equation:
1u + k2 u = 0, in D, u = fD , on Γ ,
∂u = fN , ∂n
(2.1) (2.2)
on Γ ,
(2.3)
where n is the unit normal with respect to the boundary ∂ D, and k > 0 is the wave number. In this paper, we assume that −k2 is not an eigenvalue of the Laplacian operator with the homogeneous Neumann boundary condition. Without loss of generality, we assume that the measured data fD ∈ H 1 (Γ ) and fN ∈ L2 (Γ ), and the Cauchy problem has a unique solution u in H 3/2 (D). The approximation can be split into two steps. First, we approximate the solution of the Cauchy problem by a single layer potential. Second, we approximate the single layer potential by Bessel functions. 2.2. Single layer potential approximation Define the single layer potential vϕ by
vϕ (x) :=
∂B
Φ (x, y )ϕ(y )ds(y ),
x ∈ B,
where B is a Lipschitz domain, vϕ ∈ L2 (∂ B) and Φ (x, y ) = equation. Then we have (see [16]).
(2.4) i (1) H 4 0
(k|x − y |) is the fundamental solution of the Helmholtz
Lemma 2.1. Let u ∈ H 3/2 (D) satisfy the Helmholtz equation. Let B be a bounded and simply connected domain with ∂ B ∈ C 2 , ¯ ⊂⊂ B. Then for every ε > 0, there exists a single layer potential of the form and D
vϕ (x) :=
∂B
Φ (x, y )ϕ(y )ds(y ),
x ∈ B,
(2.5)
for some ϕ ∈ L2 (∂ B), such that
∥vϕ − u∥H 3/2 (D) 6 ε,
(2.6)
∂vϕ ∂u ∥vϕ − u∥L2 (∂ D) + − 6 ε. ∂n ∂ n L2 (∂ D)
(2.7)
and
∂v
Proof. It is sufficient to prove that the set {( ∂ nϕ + ivϕ )|∂ D : ϕ ∈ L2 (∂ B)} is dense in L2 (∂ D). Let φ ∈ L2 (∂ D), satisfying
∂D
∂vϕ + ivϕ φ¯ ds(x) = 0, ∂n
for all ϕ ∈ L2 (∂ B).
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
185
Then
∂ Φ (x, y )ϕ(y )ds(y ) + i Φ (x, y )ϕ(y )ds(y ) φ(x)ds(x) ∂ D ∂ nx ∂ B ∂B ∂ Φ (x, y ) + iΦ (x, y ) φ(x)ds(x)ds(y ) = 0, = ϕ(y ) ∂ nx ∂B ∂D
for all ϕ ∈ L2 (∂ B). Define
w( ˜ y ) :=
∂D
∂ Φ (x, y ) + iΦ (x, y ) φ(x)ds(x), ∂ nx
y ∈ ∂ B,
then it is readily seen that w( ˜ y ) = 0 on ∂ B. This means that w ˜ solves the boundary value problem
1w ˜ + k2 w ˜ = 0 in R2 \ D¯ ⊃ R2 \ B¯ , w ˜ = 0 on ∂ B.
¯ And a unique continuation argument Since w ˜ satisfies Sommerfeld radiation condition, it is obvious that w( ˜ y ) = 0 in R2 \ B. 2 ¯ yields that w( ˜ y ) = 0 in R \ D. From the jump relation, we have ¯ w ˜+ −w ˜ − = φ,
∂w ˜+ ∂w ˜+ ¯ − = −iφ, ∂n ∂n
∂w ˜
on ∂ D.
∂w ˜
Since ∂ n+ + iw ˜ + = 0 on ∂ D, we have that ∂ n− + iw ˜ − = 0 on ∂ D. Therefore, w ˜ = 0 in D, and φ = 0 on ∂ D. From the denseness, we have that for every ε > 0, there exists a single layer potential vϕ of the form (2.4) for some ϕ ∈ L2 (∂ B), such that
∂vϕ ∂u 6 ε. + i v − + iu ϕ ∂n 2 ∂n L (∂ D) From the interior regularity results of the Helmholtz equation, we have
∥vϕ − u∥H 3/2 (D) 6 ε. By using the trace theorem, we obtain the estimate (2.7).
2.3. Bessel function approximation Now we approximate vϕ in D by Bessel functions. First, we introduce the Vekua operator [17–19] in the case of Helmholtz equation 1u + k2 u = 0. Let Ω ∈ C, Ω ∗ = {z ∈ C : z¯ ∈ Ω },
G(t , t ∗ , z , z ∗ ) := J0 (k (z − t )(z ∗ − t ∗ )),
∂ G(t , t ∗ , z , z ∗ ). ∂t Let v be an analytic function in Ω , then for fixed z0 ∈ Ω , we have the Vekua operator defined as z ReV[v, z0 ](z ) := Re{G(z , z¯0 , z , z¯ )v(z ) + v(t )H (z , z¯0 , z , z¯ )dt }, z = x1 + ix2 ∈ Ω . H (t , t ∗ , z , z ∗ ) := −
(2.8)
z0
And we have the following results (see [18]). Lemma 2.2 (Vekua). Let Ω be simply connected and fixed z0 ∈ Ω . Then there exists a unique analytic function v in Ω with v(z0 ) real, such that u(x1 , x2 ) = ReV[v, z0 ](z ),
z = x1 + ix2 ∈ Ω .
If Ω is bounded, it is readily seen that
∥ReV[v, z0 ](z )∥L∞ (Ω ) 6 ∥G∥∞ ∥v∥L∞ (Ω ) +
z
∥v∥L∞ (Ω ) ∥H ∥∞ d|t | 6 KV ∥v∥L∞ (Ω ) , z0
where ∥G∥∞ , ∥H ∥∞ are the supreme of these functions on Ω × Ω ∗ × Ω × Ω ∗ , and KV only depends on Ω and k.
186
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
Besides, a direct calculation gives that ReV[z n , 0](z ) =
n!2n
cos(nθ )Jn (kr ), kn n!2n ReV[iz n , 0](z ) = − n sin(nθ )Jn (kr ), k with x1 = r cos(θ ), x2 = r sin(θ ). Now, we give the approximation theorem in L∞ with Fourier–Bessel expansion. Theorem 2.1. Assume that u is the solution of the Helmholtz equation (2.1) in Bρ := {x ∈ R2 : |x| < ρ}, then there exists a Fourier–Bessel expansion of the form uN (r , θ ) =
N
cl Jl (kr )eilθ ,
(2.9)
l=−N
and we have the estimate
∥u − uN ∥L∞ (Br ) 6 C1 τ −N ,
(2.10)
∥∇(u − uN )∥L∞ (Br ) 6 C2 τ −N ,
(2.11)
where r < R < ρ , τ = R/r > 1, and C1 , C2 are constants. Proof. By Lemma 2.2, there exists a unique analytic function v ∈ Bρ , such that u = ReV[v, 0](z ). Then there exists polynomials pN of degree N, such that [15,19,18]
∥v − pN ∥L∞ (Br ) 6 c˜1 r −N ∥v∥L∞ (∂ Bρ ) ,
(2.12)
∥∇(v − pN )∥L∞ (Br ) 6 c˜2 r −N ∥v∥L∞ (∂ Bρ ) ,
(2.13)
where c˜1 and c˜2 only depends on r and ρ . Then by using the equality J−l (kr ) = (−1)l Jl (kr ), l = 1, 2, . . . , we have
∥u − uN ∥L∞ (Br )
N ilθ cl Jl (kr )e = u − ∞ l=−N L (Br ) N N (2) (1) al Jl (kr ) sin(lθ ) al Jl (kr ) cos(lθ ) + = u − ∞ l=1 l =0 L (Br ) N N (2) (1) a˜ l iz l a˜ l z l + 6 KV v − ∞ l =1 l=0 L
6 c˜1 KV τ (1)
(1)
−N
cl +(−1)l c−l
∥v∥L∞ (∂ Bρ ) 6 C1 τ (2)
cl −(−1)l c−l
, al = where a0 = c0 , al = 2 2i Similarly, we can obtain the estimate (2.11).
−N
,
(1)
, a˜ l
(Br )
=
kl (1) a , 2l l! l
(2)
a˜ l
= − 2kl l! a(l 2) , l = 1, . . . , N. l
Thus, we obtain the main result in this section. Theorem 2.2. Let u ∈ H 3/2 (D) satisfies Helmholtz equation. Then for every ε > 0, there exists a Fourier–Bessel expansion uN of the form (2.9), such that
∥u − uN ∥H 3/2 (D) 6 ε,
(2.14)
∂ u ∂ uN − ∥u − uN ∥L2 (∂ D) + 6 ε. ∂n ∂ n L2 (∂ D)
(2.15)
and
¯ ⊂⊂ Br ⊂ Bρ ⊂ B, where r < R < ρ . Proof. Choose B a bounded and simply connected domain with ∂ B ∈ C 2 , such that D From Lemma 2.1, we have that for every ε > 0, there exists a single layer potential vϕ of the form (2.4) for some ϕ ∈ L2 (∂ B), such that ∥vϕ − u∥H 3/2 (D) 6 ε.
(2.16)
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
187
By using Theorem 2.1 for vϕ in Bρ , then for ε > 0, there exists a Fourier–Bessel expansion uN of the form (2.9), such that
∥vϕ − uN ∥L∞ (Br ) 6 C1 τ −N 6 ε, ∥∇(vϕ − uN )∥L∞ (Br ) 6 C2 τ −N 6 ε, where τ = R/r > 1. Then it is readily seen that
∥vϕ − uN ∥L2 (∂ D) 6 C˜ 1 ε, and
∂vϕ − uN 6 C˜ 2 ε. ∂n 2 L (∂ D) From the interior regularity of the Helmholtz equation, we have that
∥vϕ − uN ∥H 3/2 (D) 6 C3 ε.
(2.17)
Then by (2.16) and (2.17),
∥u − uN ∥H 3/2 (D) 6 ∥u − vϕ ∥H 3/2 (D) + ∥vϕ − uN ∥H 3/2 (D) 6 C ε. This means (2.14) is valid. By using the trace theorem, the estimate (2.15) is obtained.
Now we derive the numerical approximation of uN . First, define the trace operator An : C2N +1 → L2 (Γ ) × L2 (Γ ) by
N
l=−N
An cn :=
cl Jl (kr )eilθ
N ∂ cl Jl (kr )e ∂ n l=−N
, ilθ
x ∈ Γ,
(2.18)
where cn = (c−N , c−N +1 , . . . , cN −1 , cN ) ∈ C2N +1 . Then we have the following theorem. Theorem 2.3. The operator An : C2N +1 → L2 (Γ ) × L2 (Γ ) defined by (2.18) is injective and compact. Proof. Let An cn = 0, then there exists a Fourier–Bessel expansion uN (x) =
N
∂ uN | ∂n Γ
l=−N
cl Jl (kr )eilθ , such that uN |Γ = 0 and
= 0. From the Green formula and the analyticity of the harmonic polynomial, it is readily seen that uN = 0 in R2 , and so cn = 0. Therefore, the operator An is injective. 2N +1 Let {c (k) }∞ , then there exists a convergent subsequence {c (kj ) }∞ k=1 be a bounded sequence in C j=1 , such that limj→∞ |c (kj ) − c˜ | = 0 for some c˜ ∈ C2N +1 .
√
Using the well known results of Bessel functions: |J0 (kr )| < 1, |Jl (kr )| < 1/ 2 for l = 1, 2, . . . , and J−l (kr ) = (−1)l Jl (kr ) for l = 1, 2, . . . , we have
∥J0 (kr )eilθ ∥2L2 (∂ D) 6 2π ,
∥Jl (kr )eilθ ∥2L2 (∂ D) 6 π ,
l = ±1, ±2, . . . ,
and realizing that Jl′ (r ) = Jl−1 (r ) + Jl+1 (r ), we have
2 ∂ Jl (kr )eilθ ∂n 2
L (∂ D)
2 ∂ ilθ 6 C Jl (kr )e ∂r 2
L (∂ D)
2 ∂ ilθ + Jl (kr )e ∂θ 2
L (∂ D)
6C
k2 2
+
l2
2
,
where C is a constant only depending on ∂ D. So by using the Cauchy inequality, we have
∥An c
(kj )
− An c˜ ∥
2 L2 (Γ )×L2 (Γ )
2 2 N N ∂ (kj ) (kj ) ilθ ilθ = (cl − c˜l )Jl (kr )e + (cl − c˜l )Jl (kr )e l=−N 2 ∂ n l=−N 2 L (Γ ) L (Γ ) 2 N ∂ ilθ 6 |c (kj ) − c˜ |2 ∥Jl (kr )eilθ ∥2L2 (∂ D) + ∂ n Jl (kr )e 2 L (∂ D) l=−N 2 k N ( N + 1)(2N + 1) (1 + 2N ) + . 6 |c (kj ) − c˜ |2 2π (1 + N ) + C 2
So limj→∞ ∥An c
(kj )
− An c˜ ∥L2 (Γ )×L2 (Γ ) = 0. Therefore, An is compact.
6
188
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
3. A numerical method for solving the problem Numerically, we should solve the coefficients by the equation
(An cn )(x) = b(x),
x ∈ Γ,
(3.1)
where b(x) = (f (x), g (x))T ∈ L2 (Γ ) × L2 (Γ ). But (3.1) is a severely ill-posed problem and we need a regularization method. In practice, the measured Cauchy data b(x) is often perturbed, that is, the measured data is bδ with ∥bδ − b∥ 6 δ . In this section, we introduce the truncated singular value decomposition (T-SVD) method to overcome the ill-posedness of (3.1) with perturbed measured data, i.e., we solve the equation An cn = bδ .
(3.2)
First, we introduce the singular system of An : C → L (Γ ) × L (Γ ). Let µ1 > µ2 > · · · µ2N +1 > 0 be the singular +1 +1 values of An , and (xj )2N ⊂ C2N +1 , (yj )2N ⊂ L2 (Γ ) × L2 (Γ ) be orthonormal systems with j =1 j =1 2N +1
An xj = µj yj ,
A∗n yj = µj xj ,
2
2
j = 1, . . . , 2N + 1,
where A∗n : L2 (Γ ) × L2 (Γ ) → C2N +1 is the adjoint operator of An . Then for a chosen α = α(δ) > 0, cnα(δ),δ =
1 (bδ , yj )xj µj 2
(3.3)
µj >α
is the T-SVD solution of (3.2). And we have the following theorem. Theorem 3.1. (a) There exists a function h ∈ L2 (Γ ) × L2 (Γ ), such that for every ε > 0,
|cn − A∗n h| < ε. (b) Let ∥h∥L2 (Γ )×L2 (Γ ) 6 E, then for the choice α(δ) = c δ/E for some c > 0, we have √ √ 1 δ E + ε. |cnα(δ),δ − cn | 6 √ + c
(3.4)
c
Proof. (a) Since An : C2N +1 → L2 (Γ ) × L2 (Γ ) is injective and compact, the range of A∗n is dense in C2N +1 . Therefore, for every ε > 0, there exists a function h ∈ L2 (Γ ) × L2 (Γ ), such that
|cn − A∗n h| < ε. (b) The proof can be obtained by the combination of (a) and Theorem 2.9 in book [20]. δ
δ
Numerically, the measured data b is in the discretization form (f (x1 ), f (x2 ), . . . , f δ (xm ), g δ (x1 ), g δ (x2 ), . . . , g δ (xm ))T := β δ , with Γm := {x1 , x2 , . . . , xm } ⊂ Γ , and m large enough to approximate the space L2 (Γ ) × L2 (Γ ). Denote Qm : L2 (Γ ) × L2 (Γ ) → Cm mapping bδ to β δ , then the problem (3.2) is solved in the finite space in the form Qm An cnδ = β δ ,
δ
(3.5)
δ
with |β − β| 6 δ , where β is the solution of (3.5) without noise. Now we apply the T-SVD method [21] to the discretization case (3.5). Note that Qm An : C2N +1 → Cm is a matrix, and denote it by A. We will show that a classical T-SVD method for matrix will do for (3.5) and that the solution of (3.5) will converge to that of (3.2) when m → ∞. The well known SVD of A can be written as A = U
Σ 0
V ∗ with Um×(2N +1) and V(2N +1)×(2N +1) unitary. Denote U =
(u1 , u2 , . . . , um ), V = (v1 , v2 , . . . , v2N +1 ), Σ = diag(λ1 , λ2 , . . . , λ2N +1 ), and assume that λ1 > λ2 > · · · > λ2N +1 > 0. Let 1 −1 −1 ˜ = diag(λ− α(δ) := Λ2 > λ2K for chosen K (1 6 K 6 2N + 1), and Σ 1 , λ2 , . . . , λK , 0, . . . , 0). Then the T-SVD solution of (3.5) can be written as
˜ U ∗ β δ := A˜ β δ . cnα(δ),δ = VΣ ,m
(3.6)
And a similar error estimate to Theorem 3.1 can be obtained. Theorem 3.2. Let cn,m be the solution of (3.5) with β the righthand side, and there exists z ∈ Cm with |z | 6 E, such that cn,m = A∗ z, then
δ + ΛE . |cnα(δ),δ − cn,m | 6 ,m Λ √ If we choose Λ = c δ/E for some c > 0, then √ √ 1 α(δ),δ |cn,m − cn,m | 6 √ + c δE . c
(3.7)
(3.8)
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
189
Proof.
|cnα(δ),δ − cn,m | = |A˜ β δ − A−1 β| ,m = |A˜ β δ − A˜ β + A˜ β − A−1 β| ˜ β δ − A˜ β| + |A˜ β − A−1 β|. 6 |A It is obvious that |A˜ β δ − A˜ β| 6 Λ−1 δ , and we have that
+1 +1 2N 2N 1 1 −1 ˜ (β, uj )vj = (An cn,m , uj )vj |Aβ − A β| = j=K +1 λj j=K +1 λj +1 +1 2N 2N = (cn,m , vj )vj = (Az , vj )vj j=K +1 j=K +1 +1 +1 2N 2N = (z , A∗ vj )vj = (z , λj uj )vj j=K +1 j=K +1 +1 2N λj (z , uj )vj 6 λK +1 |z | 6 ΛE . = j=K +1 So,
δ + ΛE . Λ √ And if we choose Λ = c δ/E for some c > 0, it is readily seen that √ √ 1 + c δE . |cnα(δ),δ − c | 6 √ n ,m ,m |cnα(δ),δ − cn,m | 6 ,m
c
Remark 1. It is easy to implement the numerical method. It is nothing but the traditional SVD for the matrix A, which is the discretization form of the trace operator An . Finally, we have the convergence result. Lemma 3.1. Let the assumptions of Theorems 3.1 and 3.2 hold, then the solution of (3.5) converges to that of (3.2) with the same righthand side bδ and the same regularization parameter α(δ) when m tends to infinity, i.e. lim |cnα(δ),δ − cnα(δ),δ | = 0. ,m
m→∞
Proof. It is the direct result of Remark 2.3 in [22].
Then it is readily seen that Theorem 3.3. Let the assumptions of Theorems 3.1 and 3.2 hold, then for every ε > 0, choose Λ = m sufficiently large, we have
|cnα(δ),δ − cn | 6 ,m
1
√
√ + c
c
√ δ E + ε.
√
c δ/E for some c > 0, and
(3.9)
Since the approximation of the coefficients is obtained, the approximation of the solution of the Helmholtz equation can be written as uα(δ),δ n,m (x) =
N
α(δ),δ
cl,m
Jl (kr )eilθ ,
x ∈ R2 .
l=−N
And we have the main result in the paper.
(3.10)
190
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
(a) Re(uN ) and u.
(b) Re(∂n uN ) and ∂n u. Fig. 1. The exact solution u and the numerical solution uN on Σ with k = 5 for Example 1.
Theorem 3.4. Let the assumptions of Theorems 3.1 and 3.2 hold, then for every ε > 0, choose Λ = m sufficiently large, the following estimate on boundary Σ holds,
√
c δ/E for some c > 0, and
2 ∂ uα(δ),δ ∂u n ,m 6 (1/c + c + 2)δ E − − u∥ + ∂n ∂n 2 L (Σ ) 2 k N (N + 1)(2N + 1) · 2π (1 + N ) + C (1 + 2N ) + + ε,
∥uα(δ),δ n,m
2 L2 (Σ )
2
(3.11)
6
where C only depends on ∂ D. Proof. By using the Cauchy inequality, for any x ∈ ∂ D, we have
2 ∂ uα(δ),δ ∂u n ,m − u∥ − + ∂n ∂n 2 L (∂ D) N ∂(Jl (kr )eilθ ) 2 ilθ 2 α(δ),δ 2 ∥Jl (kr )e ∥L2 (∂ D) + 6 |cn,m − cn | 2 ∂n L (∂ D) l=−N 2 k N (N + 1)(2N + 1) 6 ((1/c + c + 2)δ E + ε) 2π (1 + N ) + C (1 + 2N ) + .
∥uα(δ),δ n,m
2 L2 (∂ D)
2
So, the estimate (3.11) is obvious.
6
4. Numerical examples In this section, we present two examples to demonstrate the competitiveness of our algorithm. Example 1. Consider the case in which the exact solution to Cauchy problem is u(x) =
1 2k2 2
√
sin( 2kx2 )(ekx1 + e−kx1 )
(see [12,15]). Let D = {(x1 , x2 )|(x1 − 0.5)2 + x22 < 0.52 , x2 > 0}, Γ = {(x1 , x2 )|(x1 − 0.5) + x22 = 0.52 , x2 > 0} and Σ = ∂ D \ Γ = {(x1 , x2 )|0 6 x1 6 1, x2 = 0}. Fig. 1 shows the exact solution and the numerical solution for various levels of noise with k = 5. Table 1 shows the relative L2 error of uN and ∂ uN on Σ , respectively. From the figure and the table, we can conclude that the numerical solution is a stable approximation of the exact solution, and it converges to the exact solution when the noise decreases. Fig. 2 shows the numerical solution for various levels of wavenumber with noise level σ = 0 and σ = 5%, respectively. The result tells us that the solution for smaller k can be more stable. Example 2. Consider the unit disc D = {(x1 , x2 )|x21 + x22 < 1}. Let Γ = {x ∈ ∂ D|0 < θ (x) < Θ } and Σ = ∂ D \ Γ = {x ∈ ∂ D|Θ < θ (x) < 2π }, where θ (x) is the polar angle of x and Θ is a specified angle. In this example, we observe the effect of Θ on the numerical solution. Choose u(x) = exp(ix · d ) a plane wave, where d = (cos(π /3), sin(π /3)).
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
191
Table 1 Errors for boundary data with different noise levels for Example 1. Noise level
0 1% 3% 5%
∥uN − u∥L2 (Σ ) k=1
k=3
k=5
∥∂ uN − ∂ u∥L2 (Σ ) /∥∂ u∥L2 (Σ ) k=1 k=3
k=5
3.11e−08 3.23e−03 9.19e−03 1.37e−02
1.71e−07 3.54e−03 1.20e−02 1.40e−02
2.70e−06 7.11e−03 2.44e−02 3.37e−02
1.91e−07 9.18e−03 2.19e−02 2.85e−02
6.06e−06 4.65e−03 1.77e−02 2.73e−02
(a) Re(uN ) and u, for σ = 0.
1.09e−06 1.48e−02 4.02e−02 5.65e−02
(b) Re(uN ) and u, for σ = 5%.
Fig. 2. The exact solution u and the numerical solution uN on Σ with σ = 0 and σ = 5% for Example 1.
(a) Re(uN ) and Re(u), for σ = 5% and k = 1.
(b) Re(∂n uN ) and Re(∂n u), for σ = 5% and k = 1.
Fig. 3. The exact solution u and the numerical solution uN on Σ with k = 1, 5% noise and different values of Θ for Example 2.
Figs. 3 and 4 compare the numerical results on the whole boundary ∂ D for different values of Θ , where the wavenumber k = 1 and the noise level σ = 5% in Fig. 3, and wavenumber k = 5 and the noise level σ = 1% in Fig. 4. From these figures, it is readily seen that large Θ makes the numerical reconstruction more stable and accurate. Tables 2 and 3 show the relative L2 error of uN and ∂ uN on Σ , respectively. The results tell us that larger Θ makes the reconstruction better. 5. Conclusion In this paper, we study the Fourier–Bessel function expansion method for solving the Cauchy problem connected with the Helmholtz equation. The numerical method is based on the truncated singular value decomposition (T-SVD). Convergence
192
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
(a) Re(uN ) and Re(u), for σ = 1% and k = 5.
(b) Re(∂n uN ) and Re(∂n u), for σ = 1% and k = 5.
Fig. 4. The exact solution u and the numerical solution uN on Σ with k = 5, 1% noise and different values of Θ for Example 2. Table 2 Relative L2 errors on boundary Σ for different values of Θ and noise levels for Example 2: ∥uN − u∥L2 (Σ ) /∥u∥L2 (Σ ) . Noise level
0 1% 3% 5%
k = 0.5
k=5
Θ = π/2
Θ=π
Θ = 3π/2
Θ = π/2
Θ=π
Θ = 3π/2
7.01e−08 2.81e−02 2.07e−02 2.46e−02
1.11e−09 2.03e−03 1.62e−02 8.64e−03
9.29e−12 5.45e−04 5.12e−03 3.35e−03
4.99e−03 2.75e−01 3.46e−01 4.22e−01
2.67e−06 8.62e−02 1.23e−01 1.29e−01
1.89e−09 1.54e−03 6.97e−03 1.27e−02
Table 3 Relative L2 errors on boundary Σ for different values of Θ and noise levels for Example 2: ∥∂ uN − ∂ u∥L2 (Σ ) /∥∂ u∥L2 (Σ ) . Noise level
0 1% 3% 5%
k = 0.5
k=5
Θ = π/2
Θ=π
Θ = 3π/2
Θ = π/2
Θ=π
Θ = 3π/2
8.68e−07 1.31e−01 9.23e−02 1.04e−01
1.40e−08 1.56e−02 7.38e−02 1.19e−02
1.96e−10 4.91e−03 5.20e−02 2.80e−02
8.80e−03 2.46e−01 3.28e−01 3.59e−01
6.60e−06 1.01e−01 1.17e−01 1.20e−01
8.62e−09 4.64e−03 1.61e−02 2.35e−02
and stability are analyzed. The effectiveness of the numerical solution depends on the domain, the wavenumber, and the number of Bessel functions used. The method does not need meshing and computes fast. Acknowledgments The first author was supported by the Fundamental Research Funds for the Central Universities Grant 2014JKF03104, and the second author was supported by the National Natural Science Foundation of China (NSFC) Grant 11271159. The authors appreciate very much for the reviewers’ valuable suggestions and comments. References [1] M.R. Bai, Application of bem (boundary element method)-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries, J. Acoust. Soc. Am. 92 (1992) 533–549. [2] W. Hall, X. Mao, A boundary element investigation of irregular frequencies in electromagnetic scattering, Eng. Anal. Bound. Elem. 16 (3) (1995) 245–252. [3] B.-K. Kim, J.-G. Ih, On the reconstruction of the vibro-acoustic field over the surface enclosing an interior space using the boundary element method, J. Acoust. Soc. Am. 100 (5) (1996) 3003–3016. [4] L. Marin, L. Elliott, P. Heggs, D. Ingham, D. Lesnic, X. Wen, An alternating iterative algorithm for the Cauchy problem associated to the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 192 (5) (2003) 709–722. [5] W. Chen, Z.-J. Fu, Boundary particle method for inverse Cauchy problems of inhomogeneous Helmholtz equations, J. Mar. Sci. Technol. 17 (3) (2009) 157–163. [6] E.R.G. Alessandrini, L. Rondi, S. Vessella, The stability for the Cauchy problem for elliptic equations, Inverse Problems 25 (3) (2009) 1–47. [7] V.I.T. Hrycak, Increased stability in the continuation of solutions to the Helmholtz equation, Inverse Problems 20 (3) (2004) 697–712.
M. Liu et al. / Journal of Computational and Applied Mathematics 311 (2017) 183–193
193
[8] V. Isakov, Inverse Problems for Partial Differential Equations, in: Applied Mathematical Sciences, vol. 120, Springer-Verlag, New York, 1998, http://dx.doi.org/10.1007/978-1-4612-5338-9. [9] B. Jin, Y. Zheng, A meshless method for some inverse problems associated with the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 195 (19) (2006) 2270–2288. [10] L. Marin, D. Lesnic, The method of fundamental solutions for the Cauchy problem associated with two-dimensional Helmholtz-type equations, Comput. Struct. 83 (4) (2005) 267–278. [11] T. Wei, Y. Hon, L. Ling, Method of fundamental solutions with regularization techniques for Cauchy problems of elliptic operators, Eng. Anal. Bound. Elem. 31 (4) (2007) 373–385. [12] T. Wei, H. Qin, R. Shi, Numerical solution of an inverse 2D Cauchy problem connected with the Helmholtz equation, Inverse Problems 24 (3) (2008) 035003. [13] A. Zakeri, A. Aminataei, Q. Jannati, Application of He’s Homotopy perturbation method for Cauchy problem of ill-posed nonlinear diffusion equation, Discrete Dyn. Nat. Soc. (2010) http://dx.doi.org/10.1155/2010/780207. [14] M.A. Jafari, A. Aminataei, Exact solution for the inverse problem of finding a source term in a parabolic equation, J. Stat. Math. 2 (2) (2011) 37–39. [15] D. Zhang, F. Ma, E. Zheng, A Herglotz wave function method for the inverse Cauchy problem of the Helmholtz equation, J. Comput. Appl. Math. 237 (1) (2013) 215–222. [16] Y. Sun, D. Zhang, F. Ma, A potential function method for the Cauchy problem of elliptic operators, J. Math. Anal. Appl. 395 (1) (2012) 164–174. [17] I. Vekua, New Methods of Solution of Elliptic Equations, Gostekhizdat, Moscow, 1948. [18] T. Betcke, A gsvd formulation of a domain decomposition method forplanar eigenvalue problems, IMA J. Numer. Anal. 27 (3) (2007) 451–478. [19] J. Melenk, Operator adapted spectral element methods i: harmonic and generalized harmonic polynomials, Numer. Math. 84 (1) (1999) 35–69. [20] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, in: Applied Mathematical Sciences, vol. 120, Springer-Verlag, New York, 1996, http://dx.doi.org/10.1007/978-1-4612-5338-9. [21] M. Liu, F. Ma, Numerical method for inverse scattering in two-layered background in near-field optics, Inverse Probl. Sci. Eng. 23 (1) (2015) 130–154. http://dx.doi.org/10.1080/17415977.2014.890611. [22] C. Vogel, Optimal choice of a truncation level for the truncated svd solution of linear first kind integral equations when data are noisy, SIAM J. Numer. Anal. 23 (1) (1986) 109–117.