A Wavelet Method for the Cauchy Problem for the Helmholtz Equation

Abstract

We consider a Cauchy problem for the Helmholtz equation at a fixed frequency. The problem is severely ill posed in the sense that the solution (if it exists) does not depend continuously on the data. We present a wavelet method to stabilize the problem. Some error estimates between the exact solution and its approximation are given, and numerical tests verify the efficiency and accuracy of the proposed method.

1. Introduction

The Helmholtz equation is often used to approximate model wave propagation in inhomogeneous media. The demand for reliable numerical solutions to such type of problems is frequently encountered in geophysical and optoelectronic applications [1, 2]. In geophysical applications, for example, wave propagation simulations are used for the development of acoustic imaging techniques for gaining knowledge about geophysical structures deep within the Earth’s subsurface [3]. In optoelectronics, the determination of a radiation field surrounding a source of radiation (e.g., a light emitting diode) is also a frequently occurring problem [4]. In many engineering problems, the boundary conditions are often incomplete, either in the form of underspecified and overspecified boundary conditions on different parts of the boundary or the solution is prescribed at some internal points in the domain. These so-called Cauchy problems are inverse problems, and it is well known that they are generally ill posed in the sense of Hadamard [5]. However, the Cauchy problem suffers from the nonexistence and instability of the solution.

In this paper we consider the Cauchy problem for the Helmholtz equation in a “strip” 0<𝑥<1 as follows:Δ𝑢(𝑥,𝑦)+𝑘²𝑢𝑢(𝑥(,0𝑦,)𝑦)==0,𝑔(𝑦𝑥),∈𝑦(0∈,1ℝ),_𝑛𝑦∈ℝ^𝑛,𝑛≥1,𝑢_𝑥,(0,𝑦)=0,𝑦∈ℝ^𝑛,(1.1) where Δ=𝜕²/𝜕𝑥²+∑^𝑛𝑖=1𝜕²/𝜕𝑦^2𝑖 is an 𝑛+1 dimensional Laplace operator. We want to determine the solution 𝑢(𝑥,𝑦) for 0<𝑥≤1 from the data 𝑔(𝑦). Due to the importance of its application, this problem has been studied by many researchers, for example, DeLillo et al. [6, 7], Jin and Zheng [8], Johansson and Martin [9], and Marin et al. [10–14].

Let 𝒮 be the Schwartz space over ℝ^𝑛, and let 𝒮′ be its dual (the space of tempered distributions). Let 𝑓 denote the Fourier transform of function 𝑓(𝑦)∈𝒮 defined by𝑓(𝜉)=1(2𝜋)^𝑛/2_ℝ𝑛𝑒^{−𝑖𝜉⋅𝑦}𝑓(𝑦)𝑑𝑦,𝜉=𝜉₁,…,𝜉_𝑛,𝑦=𝑦₁,…,𝑦_𝑛,(1.2) while the Fourier transform of a tempered distribution 𝑓∈𝒮′ is defined by𝑓,𝜙=𝑓,𝜙,∀𝜙∈𝒮.(1.3) In this paper, we will consider functions depending on the variables 𝑥∈[0,1], 𝑦∈ℝ^𝑛.

For 𝑠∈ℝ, the Sobolev space 𝐻^𝑠(ℝ^𝑛) consists of all tempered distributions 𝑓(𝑦)∈𝒮′, for which 𝑓(𝜉)(1+|𝜉|²)^𝑠/2 is a function in 𝐿²(ℝ^𝑛). The norm on this space is given by‖𝑓‖_𝐻^𝑠∶=_ℝ𝑛||𝑓(𝜉)||²(1+||𝜉||²)^𝑠𝑑𝜉^1/2.(1.4)

We assume there exists a unique solution 𝑢(𝑥,𝑦) of problem (1.1), which satisfies the problem in the classical sense and 𝑔(⋅), 𝑢(𝑥,⋅)∈𝐿²(ℝ^𝑛). Applying the Fourier transform technique to problem (1.1) with respect to the variable 𝑦 yields the following problem in the frequency space:𝑢̂_𝑥𝑥(𝑥,𝜉)+𝑘²−||𝜉𝑢||²̂(0𝑢,𝜉̂()𝑥,=𝜉)𝑔̂=(𝜉0),,𝑥𝜉∈∈ℝ(0_𝑛,1),𝜉∈ℝ^𝑛,𝑛≥1,𝑢_𝑥̂,(0,𝜉)=0,𝜉∈ℝ^𝑛.(1.5) It is easy to obtain the solution of problem (1.5) (if exists) has the form𝑢̂(𝑥,𝜉)=cosh𝑥||𝜉||²−𝑘²𝑔̂(𝜉),(1.6) or equivalently, the solution of problem (1.1) has the representation𝑢(𝑥,𝑦)=1(2𝜋)^𝑛/2_ℝ𝑛𝑒^{𝑖𝜉⋅𝑦}cosh𝑥||𝜉||²−𝑘²𝑔̂(𝜉)𝑑𝜉.(1.7) Since cosh(𝑥√|𝜉|²−𝑘²) increases rapidly with exponential order as |𝜉|→∞, the Fourier transform of the exact data 𝑔(𝑦) must decay rapidly. However, in practice, the data at 𝑥=0 is often obtained on the basis of reading of physical instrument which is denoted by 𝑔_𝑚. We assume that 𝑔(⋅) and 𝑔_𝑚(⋅) satisfy‖‖𝑔(⋅)−𝑔_𝑚(⋅)‖‖_𝐻𝑟≤𝛿.(1.8) Since 𝑔_𝑚(⋅) belong to 𝐿²(ℝ^𝑛)⊂𝐻^𝑟(ℝ^𝑛) for 𝑟≤0, 𝑟 should not be positive. A small perturbation in the data 𝑔(𝑦) may cause a dramatically large error in the solution 𝑢(𝑥,𝑦) for 0<𝑥≤1. Hence problem (1.1) is severely ill posed and its numerical simulation is very difficult. It is obvious that the ill-posedness of the problem is caused by the perturbation of high frequencies.

By (1.6) we know𝑢̂(1,𝜉)=cosh||𝜉||²−𝑘²𝑔̂(𝜉).(1.9) Since the convergence rates can only be given under a priori assumptions on the exact solution [15], we will formulate such an a priori assumption in terms of the exact solution at 𝑥=1 by considering‖𝑢(1,⋅)‖_𝐻^𝑠≤𝐸.(1.10)

Meyer wavelets are special because, unlike most other wavelets, they have compact support in the frequency domain but not in the time domain (however, they decay very fast). The wavelet methods have been used to solve one-dimensional heat conduction problems [16, 17] and noncharacteristic Cauchy problem for parabolic equation in one-dimensional [18] and multidimensional [19] cases, and so forth. In this paper we propose a similar wavelet method as suggested in [19] to the problem (1.1).

The paper is organized as follows. In Section 2 we describe the Meyer wavelets and discuss the properties that make them useful for solving ill-posed problems. Some error estimates between the exact solution and its approximation as well as the choice of the regularization parameter are given in Section 3. Finally, in Section 4 numerical tests verify the efficiency and accuracy of the proposed method.

2. The Meyer Wavelets

In the present paper let Φ be Meyer’s orthonormal scaling function in 𝑛 dimensions. This function is constructed from the one-dimensional scaling functions in the following way. Let 𝜙(𝑥) and 𝜓(𝑥) be the Meyer scaling and wavelet function in one dimension defined by their Fourier transform in [20] which satisfysupp𝜙=−4supp𝜓=−83𝜋,34𝜋,3𝜋,−23𝜋∪23,83𝜋.(2.1) It can be proved (cf. [20]) that the set of functions𝜓_𝑗𝑘(𝑥)=2^𝑗/2𝜓2^𝑗𝑥−𝑘,𝑗,𝑘∈ℤ,(2.2) is an orthonormal basis of 𝐿²(ℝ). Consequently, the MRA {𝑉_𝑗}_𝑗∈ℤ of Meyer is generated by𝑉_𝑗={𝜙_𝑗𝑘,𝑘su∈ppℤ},𝜙_𝑗𝑘𝜙_𝑗𝑘∶=2^𝑗/2𝜙2^𝑗𝑥−𝑘,𝑗,𝑘∈ℤ,=𝜉;||𝜉||≤43𝜋2^𝑗.(2.3)

For the construction of an 𝑛-dimensional MRA, we take tensor products of the spaces 𝑉_𝑗 (see [21, 22]). Then the scaling function Φ is given byΦ(𝑥)=^𝑛_𝑘=1𝜙𝑥_𝑘,𝑥∈ℝ^𝑛,(2.4) and any basis function Ψ in 𝑊_𝐽 can be written in the formΨ(𝑥)=2^𝑛𝐽/2𝜓2^𝐽𝑥_𝑖−𝑘_𝑖⋅_𝑚≠𝑖𝜃_𝑚2^𝐽𝑥_𝑚−𝑘_𝑚,𝑥∈ℝ^𝑛,(2.5) where 𝑘∈ℤ^𝑛, and for any 𝑚∈{1,…,𝑛}, 𝜃_𝑚 stands for 𝜙 or 𝜓. Hence we obtain from (2.1) thatsuppΦ=−43𝜋,43𝜋^𝑛𝑓(𝜉)=0for‖𝜉‖_∞,(2.6)≤23𝜋2^𝐽,𝑓∈𝑊_𝐽,𝐽∈ℕ.(2.7) The orthogonal projection on the space 𝑉_𝐽 is defined by𝑃_𝐽𝑓∶=_{𝑘∈ℤ𝑛}𝑓,Φ_𝐽,𝑘Φ_𝐽,𝑘,(2.8) while 𝑄_𝐽𝑓 denotes the orthogonal projection of a function 𝑓 on the wavelet space 𝑊_𝐽 with 𝑉_𝐽+1=𝑉_𝐽⊕𝑊_𝐽. (In many contexts one will find more than one detailed space 𝑊_𝐽, that is, 𝑉_𝐽+1=𝑉_𝐽⊕𝑊_1,𝐽⊕𝑊_2,𝐽⊕⋯. Here, the space 𝑊_𝐽 is simply defined as the orthogonal complement of 𝑉_𝐽 in 𝑉_𝐽+1).

LetΩ_𝐽∶=2^𝐽−23𝜋,32𝜋^𝑛.(2.9) Setting Γ_𝐽∶=ℝ^𝑛⧵Ω_𝐽, together with (2.6), it follows for 𝐽∈ℕ that𝐼−𝑃_𝐽𝑃_𝐽𝑓(𝜉)=0for𝜉∈Γ_𝐽+1,𝑓̂(𝜉)=𝑄_𝐽𝑓(𝜉)for𝜉∈Ω_𝐽+1.(2.10) We introduce the operator 𝑀_𝐽 which is defined by the equation𝑀_𝐽𝑓∶=1−𝜒_𝐽𝑓,𝐽∈ℕ,(2.11) where 𝜒_𝐽 denotes the characteristic function of the cube Ω_𝐽. From (2.7) it follows that any basis function Ψ in 𝑊_𝑗, 𝑗≥𝐽, satisfiesΨ(𝜉)=0,𝜉∈Ω_𝐽,(2.12) and we obtain(𝑓,Ψ)=𝑓,Ψ=1−𝜒_𝐽𝑓,Ψ=𝑀_𝐽𝑓,Ψ.(2.13) And it follows for 𝐽∈ℕ that𝐼−𝑃𝑄_𝐽𝐽=𝑄_𝐽𝑀_𝐽,=𝐼−𝑃_𝐽𝑀_𝐽.(2.14)

3. Wavelet Regularization and Error Estimates

We list the following two lemmas given in [19, 23] which are useful to our proof.

Lemma 3.1 (see [19, 23]). Let {𝑉_𝐽}_𝐽∈ℤ be an m-regular MRA, and let 𝑟,𝑠∈ℝ be such that −𝑚<𝑟<𝑠<𝑚. Then for each function 𝑓∈𝐻^𝑠(ℝ^𝑛) and 𝐽∈ℕ, the following inequality holds: ‖‖𝑓−𝑃_𝐽𝑓‖‖_𝐻𝑟≤𝐶₁2^{−𝐽(𝑠−𝑟)}‖𝑓‖_𝐻^𝑟.(3.1)

Lemma 3.2 (see [19]). Let {𝑉_𝐽}_𝐽∈ℤ be Meyer’s (tensor-) MRA, and suppose 𝐽∈ℕ, 𝑟∈ℝ. Then for all 𝑓∈𝑉_𝐽, one has ‖‖‖𝜕^𝑙𝜕𝑥^𝑙𝑖𝑓‖‖‖_𝐻𝑟≤𝐶₂2^{(𝐽−1)𝑙}‖𝑓‖_𝐻^𝑟,𝑖=1,…,𝑛,𝑙∈ℕ.(3.2)

Define an operator 𝑇_𝑥∶𝑔(𝑦)↦𝑢(𝑥,𝑦) by (1.6), that is,𝑇_𝑥𝑔=𝑢(𝑥,𝑦),0<𝑥≤1,(3.3) or equivalently,𝑇_𝑥𝑔(𝜉)=cosh𝑥||𝜉||²−𝑘²𝑔̂(𝜉),0<𝑥≤1.(3.4) Then we have

Theorem 3.3. Let {𝑉_𝐽}_𝐽∈ℤ be Meyer’s MRA and suppose 𝑟∈ℝ and 𝐽∈ℕ which satisfies 2^𝐽>𝑘, 0≤𝑥≤1. Then for all 𝑓∈𝑉_𝐽, one has ‖‖𝑇_𝑥𝑓‖‖_𝐻𝑟≤𝐶₅𝑒^{𝑥√22(𝐽−1)−𝑘2}+1‖𝑓‖_𝐻^𝑟.(3.5)

Proof. For 𝑓∈𝑉_𝐽, by definition (1.4) and formula (3.4), from Lemma 3.2, we have ‖‖𝑇_𝑥𝑓‖‖_𝐻𝑟=_ℝ𝑛||||cosh(𝑥||𝜉||²−𝑘²)𝑓||||²1+||𝜉||²^𝑟𝑑𝜉^1/2≤⎛⎜⎜⎜⎜⎜⎝_|𝜉|>𝑘||||||||cosh𝑥||𝜉||²−𝑘²cosh𝑥||𝜉||cosh𝑥||𝜉||𝑓||||||||²1+||𝜉||²^𝑟𝑑𝜉⎞⎟⎟⎟⎟⎟⎠^1/2+_{|𝜉|≤𝑘}||||cos𝑥𝑘²−||𝜉||²𝑓||||²1+||𝜉||²^𝑟𝑑𝜉^1/2≤sup_{𝜉∈Ω𝐽+1}||||||||cosh𝑥||𝜉||²−𝑘²cosh𝑥||𝜉||||||||||_|𝜉|>𝑘||cosh𝑥||𝜉||𝑓||²1+||𝜉||²^𝑟𝑑𝜉^1/2+‖𝑓‖_𝐻^𝑟≤2sup_{𝜉∈Ω𝐽+1}𝑒^{√𝑥(|𝜉|2−𝑘2−|𝜉|)}⎛⎜⎜⎝_|𝜉|>𝑘|||||_∞_𝑙=0𝑥^2𝑙(2𝑙)!|𝜉|^2𝑙𝑓|||||²1+||𝜉||²^𝑟𝑑𝜉⎞⎟⎟⎠^1/2+‖𝑓‖_𝐻^𝑟≤2𝐶₃𝑒^{√𝑥(22(𝐽−1)−𝑘2−2𝐽−1)∞}_𝑙=0𝑥^2𝑙(2𝑙)!‖‖(Δ^𝑦)^𝑙𝑓‖‖_𝐻𝑟+‖𝑓‖_𝐻^𝑟≤𝐶₄𝑒^{√𝑥(22(𝐽−1)−𝑘2−2𝐽−1)∞}_𝑙=0𝑥^2𝑙(2𝑙)!⋅𝑛2^{2(𝐽−1)𝑙}+1‖𝑓‖_𝐻^𝑟≤𝐶₅𝑒^{√𝑥(22(𝐽−1)−𝑘2−2𝐽−1)}cosh𝑥2^𝐽−1+1‖𝑓‖_𝐻^𝑟≤𝐶₅𝑒^{𝑥√22(𝐽−1)−𝑘2}+1‖𝑓‖_𝐻^𝑟.(3.6)

Since the Cauchy data are given inexactly by 𝑔_𝑚, we need a stable algorithm to approximate the solution of (1.1). Our method is as follows. Consider the operator𝑇_𝑥,𝐽∶=𝑃_𝐽𝑇_𝑥𝑃_𝐽,(3.7) and show that it approximates 𝑇_𝑥 in a stable way for an appropriate choice for 𝐽∈ℕ depending on 𝛿 and 𝐸. By the triangle inequality we know‖‖𝑇_𝑥𝑔−𝑇_𝑥,𝐽𝑔_𝑚‖‖_𝐻𝑟≤‖‖(𝑇_𝑥−𝑇_𝑥,𝐽)𝑔‖‖_𝐻𝑟+‖‖𝑇_𝑥,𝐽𝑔−𝑔_𝑚‖‖_𝐻^𝑟.(3.8)

From (1.8) and Theorem 3.3, the second term on the right-hand side of (3.8) satisfies‖‖𝑇_𝑥,𝐽𝑔−𝑔_𝑚‖‖_𝐻^𝑟=‖‖𝑃_𝐽𝑇_𝑥𝑃_𝐽𝑔−𝑔_𝑚‖‖_𝐻^𝑟≤‖‖𝑇_𝑥𝑃_𝐽(𝑔−𝑔_𝑚)‖‖_𝐻𝑟≤𝐶₅𝑒^{𝑥√22(𝐽−1)−𝑘2}+1𝛿.(3.9)

For the first one we have‖‖𝑇_𝑥−𝑇_𝑥,𝐽𝑔‖‖_𝐻^𝑟≤‖‖𝐼−𝑃_𝐽𝑇_𝑥𝑔‖‖_𝐻𝑟+‖‖𝑃_𝐽𝑇_𝑥𝐼−𝑃_𝐽𝑔‖‖_𝐻^𝑟.(3.10)

By Lemma 3.1, (1.9), (1.10), (2.14), and (3.4), we get‖‖𝐼−𝑃_𝐽𝑇_𝑥𝑔‖‖_𝐻𝑟=‖‖𝐼−𝑃_𝐽𝑀_𝐽𝑇_𝑥𝑔‖‖_𝐻𝑟≤𝐶₁2^{−𝐽(𝑠−𝑟)}‖‖𝑀_𝐽𝑇_𝑥𝑔‖‖_𝐻𝑠=𝐶₁2^{−𝐽(𝑠−𝑟)}⎛⎜⎜⎜⎜⎜⎝_Γ𝐽||||||||cosh𝑥||𝜉||²−𝑘²cosh||𝜉||²−𝑘²𝑢(1,⋅)||||||||²1+||𝜉||²^𝑠⎞⎟⎟⎟⎟⎟⎠^1/2≤𝐶₁2^{−𝐽(𝑠−𝑟)}sup_{𝜉∈Γ𝐽}||||||||cosh𝑥||𝜉||²−𝑘²cosh||𝜉||²−𝑘²||||||||⋅‖𝑢(1,⋅)‖^𝐻𝑠≤2𝐶₁23𝜋2^𝐽²−𝑘²^{−(𝑠−𝑟)/2}≤2𝐶₁𝑒^{−(1−𝑥)((2/3)𝜋2𝐽)2−𝑘2}𝐸2^2𝐽−𝑘²^{−(𝑠−𝑟)/2}𝑒^{√−(1−𝑥)22𝐽−𝑘2}𝐸.(3.11) On the other hand, due to (2.10), we know‖‖𝑃_𝐽𝑇_𝑥𝐼−𝑃_𝐽𝑔‖‖_𝐻^𝑟≤‖‖𝑇_𝑥𝐼−𝑃_𝐽𝑔‖‖_𝐻^𝑟≤_Ω𝐽+1||||cosh𝑥||𝜉||²−𝑘²𝑄^𝐽𝑔̂(𝜉)||||²1+||𝜉||²^𝑟𝑑𝜉^1/2+_Γ𝐽+1||||cosh𝑥||𝜉||²−𝑘²𝑔̂(𝜉)||||²1+||𝜉||²^𝑟𝑑𝜉^1/2=∶𝐼₁+𝐼₂.(3.12) We estimate the two parts at the right-hand side of (3.12) separately. For 𝐼₂ we have𝐼₂=_Γ𝐽+1||||cosh𝑥||𝜉||²−𝑘²𝑔̂(𝜉)||||²1+||𝜉||²^𝑟𝑑𝜉^1/2=⎛⎜⎜⎜⎜⎜⎝_Γ𝐽+1||||||||cosh𝑥||𝜉||²−𝑘²cosh||𝜉||²−𝑘²𝑢(1,⋅)||||||||²1+||𝜉||²^𝑟𝑑𝜉⎞⎟⎟⎟⎟⎟⎠^1/2≤sup_{𝜉∈Γ𝐽+1}||||||||cosh𝑥||𝜉||²−𝑘²cosh||𝜉||²−𝑘²1+||𝜉1||²^{(𝑠−𝑟)/2}||||||||×_Γ𝐽+1||||𝑢(1,⋅)1+||𝜉||₂^𝑠/2||||²𝑑𝜉^1/2≤2⋅43𝜋2^𝐽^{−(𝑠−𝑟)}𝑒^{−(1−𝑥)𝑛((4/3)𝜋2𝐽)2−𝑘2}⋅‖𝑢(1,⋅)‖_𝐻^𝑠≤2⋅(2^2𝐽−𝑘²)^{−(𝑠−𝑟)/2}𝑒^{√−(1−𝑥)22𝐽−𝑘2}𝐸.(3.13) Now we turn to 𝐼₁. There holds𝐼₁=_Ω𝐽+1||||cosh𝑥||𝜉||²−𝑘²𝑄^𝐽𝑔̂(𝜉)||||²1+||𝜉||²^𝑟𝑑𝜉^1/2≤‖‖𝑇_𝑥𝑄_𝐽𝑔‖‖_𝐻𝑟≤𝐶₅𝑒^{𝑥√22𝐽−𝑘2}+1‖‖𝑄_𝐽𝑔‖‖_𝐻𝑟,(3.14) since 𝑄_𝐽𝑔∈𝑉_𝐽+1. Furthermore, from (2.14), it follows that‖‖𝑄_𝐽𝑔‖‖_𝐻𝑟=‖‖𝑄_𝐽𝑀_𝐽𝑔‖‖_𝐻𝑟≤‖‖𝑀_𝐽𝑔‖‖_𝐻𝑟=_Γ𝐽||𝑔̂(𝜉)||²1+||𝜉||²^𝑟𝑑𝜉^1/2=⎛⎜⎜⎜⎜⎜⎝_Γ𝐽||||||||𝑢̂(1,𝜉)cosh||𝜉||²−𝑘²||||||||²1+||𝜉||²^𝑟𝑑𝜉⎞⎟⎟⎟⎟⎟⎠^1/2≤2⋅2^{−𝐽(𝑠−𝑟)}𝑒^{−𝑛((2/3)𝜋2𝐽)2−𝑘2}‖𝑢(1,⋅)‖_𝐻^𝑠≤2⋅(2^2𝐽−𝑘²)^{−(𝑠−𝑟)/2}𝑒^{−√22𝐽−𝑘2}𝐸.(3.15) Therefore,‖‖𝑃_𝐽𝑇_𝑥(𝐼−𝑃_𝐽)𝑔‖‖_𝐻𝑟≤+2212+_2𝐽𝐶₅2^2𝐽−𝑘²^{−(𝑠−𝑟)/2}𝑒^{√−(1−𝑥)22𝐽−𝑘2}𝐸−𝑘²^{−(𝑠−𝑟)/2}𝑒^{−√22𝐽−𝑘2}𝐸.(3.16) Combining (3.11) and (3.16) with (3.10), we have‖‖(𝑇_𝑥−𝑇_𝑥,𝐽)𝑔‖‖_𝐻𝑟≤+22𝐶2₁+_2𝐽1+𝐶₅2^2𝐽−𝑘²^{−(𝑠−𝑟)/2}𝑒^{√−(1−𝑥)22𝐽−𝑘2}𝐸−𝑘²^{−(𝑠−𝑟)/2}𝑒^{−√22𝐽−𝑘2}𝐸.(3.17) Then from (3.9) and (3.17) we finally arrive at‖‖𝑇_𝑥𝑔−𝑇_𝑥,𝐽𝑔_𝑚‖‖_𝐻𝑟≤+2𝐶₅𝑒𝐶^{𝑥√22(𝐽−1)−𝑘2}₁+1𝛿+22^2𝐽−𝑘²^{−(𝑠−𝑟)/2}𝑒^{−√22𝐽−𝑘2}𝐸+1+𝐶₅2^2𝐽−𝑘²^{−(𝑠−𝑟)/2}𝑒^{√−(1−𝑥)22𝐽−𝑘2}𝐸.(3.18)

In order to show some stability estimates of the Hölder type for our method using (3.18), we use the following lemma which appeared in [24] for choosing a proper regularization parameter 𝐽.

Lemma 3.4. Let the function 𝑓(𝜆)∶[0,𝑎]→ℝ be given by 𝑓(𝜆)=𝜆^𝑏𝑑ln1𝜆^−𝑐(3.19) with a constant 𝑐∈ℝ and positive constants 𝑎<1, 𝑏, and 𝑑. Then for the inverse function 𝑓⁻¹(𝜆), one has 𝑓⁻¹(𝜆)=𝜆^1/𝑏𝑑𝑏ln𝜆1^𝑐/𝑏(1+𝑜(1))for𝜆⟶0.(3.20)

Based on this lemma, we can choose the regularization parameter 𝐽 by minimizing the right-hand side of (3.18).

Denote𝑒^{−√22𝐽−𝑘2}=𝜆∈(0,1),(3.21) and let 𝐶=𝐶₅/2(1+𝐶₁+𝐶₅) and𝐶𝜆^−𝑥𝛿=𝜆^1−𝑥ln1𝜆^{−(𝑠−𝑟)}𝐸,(3.22) that is,𝐶𝛿𝐸=𝜆ln𝜆1^{−(𝑠−𝑟)}.(3.23) Then by Lemma 3.4 we obtain that𝜆=𝐶𝛿𝐸ln𝐶𝛿1/𝐸^𝑠−𝑟(1+𝑜(1))for𝐶𝛿=𝐶𝛿𝐸⟶0𝐸ln𝐶𝐸𝛿^𝑠−𝑟(1+𝑜(1))for𝛿⟶0.(3.24) Taking the principal part of 𝜆, we get𝐽^∗=12log²ln²𝐶𝐸𝛿ln𝐶𝐸𝛿^{−(𝑠−𝑟)}+𝑘²,(3.25) due to (3.21). Now, summarizing above inference process, we obtain the main result of the present paper.

Theorem 3.5. For 𝑠≥𝑟, suppose that conditions (1.8) and (1.10) hold. If one takes 𝐽^∗∗=𝐽^∗,(3.26) where 𝐽^∗ was defined in (3.25), [𝑎] with square bracket denotes the largest integer less than or equal to 𝑎∈𝑅. Then there holds the following stability estimate: ‖‖𝑇_𝑥𝑔−𝑇_𝑥,𝐽𝑔_𝑚‖‖≤2𝐶₁+1+𝐶₅𝐸^𝑥𝐶₅𝛿^1−𝑥ln𝐸𝐶𝛿^{−(𝑠−𝑟)𝑥}×1+ln𝐸/𝐶𝛿+1+2𝐶ln𝐸/𝐶𝛿+lnl(nln𝐸𝐸/𝐶/𝐶𝛿𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟=2𝐶₁ln𝐸/𝐶𝛿+ln(ln𝐸/𝐶𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟𝛿+1+𝐶₅𝐸^𝑥𝐶₅𝛿^1−𝑥ln𝐸𝐶𝛿^{−(𝑠−𝑟)𝑥}(1+𝑜(1)),(3.27) for 𝛿→0.

Remark 3.6. In general, the a priori bound 𝐸 and the coefficients 𝐶₁-𝐶₅ and 𝐶 are not exactly known in practice. In this case, with 𝐽^⋆=12log²ln²𝛿1ln𝛿1^{−(𝑠−𝑟)}+𝑘²,(3.28) it holds that ‖‖𝑇_𝑥𝑔−𝑇_𝑥,𝐽𝑔_𝑚‖‖≤𝛿^1−𝑥ln1𝛿^{−(𝑠−𝑟)𝑥}𝐶⁵+2𝐶₁+1+𝐶₅𝐸ln1/𝛿+𝐶₅ln1/𝛿+ln(ln1/𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟+2𝐸ln1/𝛿=𝛿_1−𝑥ln1/𝛿+ln(ln1/𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟𝛿ln𝛿1^{−(𝑠−𝑟)𝑥}(1+𝑜(1)),(3.29) for 𝛿→0.

Remark 3.7. The proposed wavelet method can also be used to solve the following Cauchy problem for the modified Helmholtz equation (i.e., the Yukawa equation [25]) Δ𝑣(𝑥,𝑦)+𝑘²𝑣𝑣(𝑥(0,,𝑦𝑦))==0𝑔,(𝑦𝑥),∈𝑦(0∈,1ℝ),_𝑛𝑦∈ℝ^𝑛,𝑛≥1,𝑣_𝑥,(0,𝑦)=0,𝑦∈ℝ^𝑛,(3.30) where Δ=𝜕²/𝜕𝑥²+∑^𝑛𝑖=1𝜕²/Δ𝑦^2𝑖 is the same as in (1.1).
It is easy to know that the exact solution of problem (3.30) is 𝑣(𝑥,𝑦)=1(2𝜋)^𝑛/2_ℝ𝑛𝑒^{𝑖𝜉⋅𝑦}cosh𝑥||𝜉||²+𝑘²𝑔̂(𝜉)𝑑𝜉.(3.31) Define an operator 𝑇_𝑥∶𝑔(𝑦)↦𝑣(𝑥,𝑦) such that 𝑇_𝑥𝑔(𝜉)=cosh𝑥||𝜉||²+𝑘²𝑔̂(𝜉),0<𝑥≤1,(3.32) and the approximate solution is 𝑣^𝛿𝐽=𝑇_𝑥,𝐽𝑔_𝑚,(3.33) where 𝑔 and 𝑔_𝑚 satisfy (1.8), and 𝑇_𝑥,𝐽=𝑃_𝐽𝑇_𝑥𝑃_𝐽. If we select the regularization parameter 𝐽^†=12log²ln²𝛿1ln𝛿1^{−(𝑠−𝑟)}−𝑘²,(3.34) then there holds ‖‖𝑣(𝑥,⋅)−𝑣^𝛿𝐽(𝑥,⋅)‖‖≤𝛿^1−𝑥ln𝛿1^{−(𝑠−𝑟)𝑥}×𝐶₅+2𝐶₁+1+𝐶₅𝐸ln1/𝛿+𝐶₅ln1/𝛿+ln(ln1/𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟+2𝐸ln1/𝛿=𝛿_1−𝑥ln1/𝛿+ln(ln1/𝛿)^{−(𝑠−𝑟)}^𝑠−𝑟𝛿ln𝛿1^{−(𝑠−𝑟)𝑥}(1+(𝑜(1))),(3.35) for 𝛿→0.

4. Numerical Aspect

4.1. Numerical Implementation

We want to discuss some numerical aspects of the proposed method in this section.

We consider the case when 𝑛=2. Supposing that the sequence {𝑔(𝑦_1,𝑖,𝑦_2,𝑗)}^{𝑁𝑖,𝑗=1} represents samples from the function 𝑔(𝑦₁,𝑦₂) on an equidistant grid in the square [𝑎,𝑏]², and 𝑁 is even, then we add a random uniformly distributed perturbation to each data and obtain the perturbation data𝑔_𝑚=𝑔+𝜇randn(size(𝑔)).(4.1) Then the total noise 𝛿 can be measured in the sense of root mean square error according to𝛿∶=‖‖𝑔_𝑚−𝑔‖‖_𝑙2=⎷1𝑁^2𝑁_𝑁𝑖=1_𝑗=1𝑔_𝑚𝑦_1,𝑖,𝑦_2,𝑗−𝑔𝑦_1,𝑖,𝑦_2,𝑗²,(4.2) where ”randn(⋅)” is a normally distributed random variable with zero mean and unit standard deviation and 𝜖 dictates the level of noise. ”randn(size(𝑔))” returns an array of random entries that is the same size as 𝑔.

For the function 𝑔_𝑚(𝑦₁,𝑦₂), we have𝑢^𝛿𝐽(𝑥,𝑦)=𝑇_𝑥,𝐽𝑔_𝑚=𝑃_𝐽𝑇_𝑥𝑃_𝐽𝑔_𝑚.(4.3) Hence, by using it with 𝐽^⋆ being given in (3.28), we can obtain the approximate solution.

We will use DMT as a short form of the “discrete Meyer (wavelet) transform." Algorithms for discretely implementing the Meyer wavelet transform are described in [21]. These algorithms are based on the fast Fourier transform (FFT), and computing the DMT of a vector in ℝ requires 𝒪(𝑁log²²𝑁) operations.

4.2. Numerical Tests

In this section some numerical tests are presented to demonstrate the usefulness of the approach. The tests were performed using Matlab and the wavelet package WaveLab 850, which was downloaded from http://www-stat.stanford.edu/~wavelab/. Throughout this section, we set 𝜇=10⁻³, 𝑎=−5, 𝑏=5, and 𝑁=2⁶.

Example 4.1. Take 𝑛=2 and 𝑔(𝑦)=𝑒^−𝑦2∈𝒮(ℝ²), where 𝑦=(𝑦₁,𝑦₂) and 𝒮(ℝ²) denotes the Schwartz function space.
Since 𝑔̂(𝜉)∈𝒮(ℝ²), 𝜉=(𝜉₁,𝜉₂) decays rapidly, and the formula (1.7) can be used to calculate 𝑢(𝑥,𝑦) with exact data directly, that is, 𝑢(𝑥,𝑦)=12𝜋_ℝ²𝑒^{𝑖(𝜉1𝑦1+𝜉2𝑦2)}cosh𝑥𝜉²¹+𝜉²²−𝑘²𝑔̂𝜉¹,𝜉₂𝑑𝜉₁𝑑𝜉₂.(4.4)

In Figure 1 we give the exact solution at 𝑥=1, that is, 𝑢(1,𝑦₁,𝑦₂), and the reconstructed solution 𝑢^𝛿(1,𝑦₁,𝑦₂) from the noisy data 𝑔_𝑚(𝑦₁,𝑦₂) without regularization. We see that 𝑢^𝛿 does not approximate the solution and some regularization procedure is necessary.

(a)

(b)

Figure 1 

(a) Exact solution 𝑢(1,⋅); (b) unregularized solution reconstructed from 𝑔𝑚 for 𝑥=1.

Letting 𝑘=1, the regularized solutions and the corresponding errors 𝑢−𝑢^𝛿𝐽_⋆ defined by the regularization parameter 𝐽^⋆=3,4,5 are illustrated in Figure 2. We can see that in 𝑉₃ the approximation is very poor since the frequencies are cut off excessively by the projection 𝑃₃. If 𝐽^⋆ is taken to be too large, the noise in the function 𝑔_𝑚 is not damped enough by 𝑃_𝐽^⋆, and thus the high frequencies of 𝑔̂_𝑚 are so extremely magnified that they destroy the approximated solution. The approximation parameter 𝐽^⋆=4 seems to be the optimal choice for this example.

Figure 2 

(a) The regularized solution at 𝑥=1; (b) difference between the regularization and the exact solution; (1), (2), (3) correspond to 𝐽⋆=3,4,5, respectively.

In Figure 3 we display the exact solution, its approximation, and corresponding errors for 𝑘=5 and 100, respectively. We see that the proposed method is useful for different wave number 𝑘.

Figure 3 

(a) and (b) correspond to 𝑘=5 and 𝑘=100, respectively; (1), (2), (3) correspond to the exact solution, the regularized solution and the difference between the regularization and the exact solution, respectively.

Figure 4 shows that the proposed method for the Cauchy problem for the modified Helmholtz equation is also effective.

(a)

(b)

(c)

(d)

Figure 4 

(a) Exact solution 𝑣(1,⋅); (b) unregularized solution reconstructed from 𝑔𝑚 for 𝑥=1; (c) regularized solution reconstructed from 𝑔𝑚 for 𝑥=1 and 𝐽†=4; (d) the difference between the regularization and the exact solution.

Acknowledgments

The authors would like to thank the WaveLab Team at the Stanford University for the help of their wavelet package Wavelab850. The work described in this paper was supported by the Fundamental Research Funds for the Central Universities of China (Project no. ZYGX2009J099) and the National Natural Science Foundation of China (Project no. 11171136).