A Wavelet Method for the Cauchy Problem for the Helmholtz Equation

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


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.

ISRN Applied Mathematics
In this paper we consider the Cauchy problem for the Helmholtz equation in a "strip" 0 < x < 1 as follows: Δu x, y k 2 u x, y 0, x ∈ 0, 1 , y ∈ R n , n ≥ 1, u 0, y g y , y ∈ R n , u x 0, y 0, y ∈ R n , where Δ ∂ 2 /∂x 2 n i 1 ∂ 2 /∂y 2 i is an n 1 dimensional Laplace operator.We want to determine the solution u x, y for 0 < x ≤ 1 from the data g y .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 S be the Schwartz space over R n , and let S be its dual the space of tempered distributions .Let f denote the Fourier transform of function f y ∈ S defined by In this paper, we will consider functions depending on the variables x ∈ 0, 1 , y ∈ R n .For s ∈ R, the Sobolev space H s R n consists of all tempered distributions f y ∈ S , for which f ξ 1 |ξ| 2 s/2 is a function in L 2 R n .The norm on this space is given by 1.4 We assume there exists a unique solution u x, y of problem 1.1 , which satisfies the problem in the classical sense and g • , u x, • ∈ L 2 R n .Applying the Fourier transform technique to problem 1.1 with respect to the variable y yields the following problem in the frequency space:

1.5
It is easy to obtain the solution of problem 1.5 if exists has the form or equivalently, the solution of problem 1.1 has the representation Since cosh x |ξ| 2 − k 2 increases rapidly with exponential order as |ξ| → ∞, the Fourier transform of the exact data g y must decay rapidly.However, in practice, the data at x 0 is often obtained on the basis of reading of physical instrument which is denoted by g m .We assume that g • and g m • satisfy Since g m • belong to L 2 R n ⊂ H r R n for r ≤ 0, r should not be positive.A small perturbation in the data g y may cause a dramatically large error in the solution u x, y for 0 < x ≤ 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 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 x 1 by considering 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.

The Meyer Wavelets
In the present paper let Φ be Meyer's orthonormal scaling function in n dimensions.This function is constructed from the one-dimensional scaling functions in the following way.Let φ x and ψ x be the Meyer scaling and wavelet function in one dimension defined by their Fourier transform in 20 which satisfy

2.1
It can be proved cf.20 that the set of functions is an orthonormal basis of L 2 R .Consequently, the MRA {V j } j∈Z of Meyer is generated by

2.3
For the construction of an n-dimensional MRA, we take tensor products of the spaces V j see 21, 22 .Then the scaling function Φ is given by and any basis function Ψ in W J can be written in the form where k ∈ Z n , and for any m ∈ {1, . . ., n}, θ m stands for φ or ψ.Hence we obtain from 2.1 that The orthogonal projection on the space V J is defined by while Q J f denotes the orthogonal projection of a function f on the wavelet space W J with V J 1 V J ⊕ W J .In many contexts one will find more than one detailed space W J , that is, Here, the space W J is simply defined as the orthogonal complement of V J in V J 1 . Let

2.9
Setting Γ J : R n \ Ω J , together with 2.6 , it follows for J ∈ N that

2.10
We introduce the operator M J which is defined by the equation where χ J denotes the characteristic function of the cube Ω J .From 2.7 it follows that any basis function Ψ in W j , j ≥ J, satisfies and we obtain And it follows for J ∈ N that 2.14

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 {V J } J∈Z be an m-regular MRA, and let r, s ∈ R be such that −m < r < s < m.Then for each function f ∈ H s R n and J ∈ N, the following inequality holds: Define an operator T x : g y → u x, y by 1.6 , that is, or equivalently, Then we have Theorem 3.3.Let {V J } J∈Z be Meyer's MRA and suppose r ∈ R and J ∈ N which satisfies Then for all f ∈ V J , one has 1 f H r .

ISRN Applied Mathematics 7
Proof.For f ∈ V J , by definition 1.4 and formula 3.4 , from Lemma 3.2, we have 1 f H r .

3.6
Since the Cauchy data are given inexactly by g m , we need a stable algorithm to approximate the solution of 1.1 .Our method is as follows.Consider the operator T x,J : P J T x P J , 3.7 and show that it approximates T x in a stable way for an appropriate choice for J ∈ N depending on δ and E. By the triangle inequality we know From 1.8 and Theorem 3.3, the second term on the right-hand side of 3.8 satisfies

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

3.11
On the other hand, due to 2.10 , we know

3.12
We estimate the two parts at the right-hand side of 3.12 separately.For I 2 we have

3.13
Now we turn to I 1 .There holds 3.14 since Q J g ∈ V J 1 .Furthermore, from 2.14 , it follows that

3.16
Combining 3.11 and 3.16 with 3.10 , we have 3.17 Then from 3.9 and 3.17 we finally arrive at

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 J.
Based on this lemma, we can choose the regularization parameter J by minimizing the right-hand side of 3. 18 .Denote

3.23
Then by Lemma 3.4 we obtain that

3.24
Taking the principal part of λ, we get where J * was defined in 3.25 , a with square bracket denotes the largest integer less than or equal to a ∈ R. Then there holds the following stability estimate:

3.27
for δ → 0. Remark 3.6.In general, the a priori bound E and the coefficients C 1 -C 5 and C are not exactly known in practice.In this case, with
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

3.30
where Δ ∂ 2 /∂x 2 n i 1 ∂ 2 /Δy 2 i is the same as in 1.1 .It is easy to know that the exact solution of problem 3.30 is v x, y 1 2π n/2 R n e iξ•y cosh x |ξ| 2 k 2 g ξ dξ.

3.31
Define an operator T x : g y → v x, y such that and the approximate solution is where g and g m satisfy 1.8 , and T x,J P J T x P J .If we select the regularization parameter

Numerical Implementation
We want to discuss some numerical aspects of the proposed method in this section.We consider the case when n 2. Supposing that the sequence {g y 1,i , y 2,j } N i,j 1 represents samples from the function g y 1 , y 2 on an equidistant grid in the square a, b 2 , and N is even, then we add a random uniformly distributed perturbation to each data and obtain the perturbation data g m g μ randn size g .

4.1
Then the total noise δ can be measured in the sense of root mean square error according to where "randn • " is a normally distributed random variable with zero mean and unit standard deviation and dictates the level of noise."randn size g " returns an array of random entries that is the same size as g.
For the function g m y 1 , y 2 , we have u δ J x, y T x,J g m P J T x P J g m .

4.3
Hence, by using it with J 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 R requires O Nlog 2 2 N operations.

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  Since g ξ ∈ S R 2 , ξ ξ 1 , ξ 2 decays rapidly, and the formula 1.7 can be used to calculate u x, y with exact data directly, that is, u x, y 1 2π R 2 e i ξ 1 y 1 ξ 2 y 2 cosh x ξ 2 1 ξ 2 2 − k 2 g ξ 1 , ξ 2 dξ 1 dξ 2 .

4.4
In Figure 1 we give the exact solution at x 1, that is, u 1, y 1 , y 2 , and the reconstructed solution u δ 1, y 1 , y 2 from the noisy data g m y 1 , y 2 without regularization.We see that u δ does not approximate the solution and some regularization procedure is necessary.
Letting k 1, the regularized solutions and the corresponding errors u−u δ J defined by the regularization parameter J 3, 4, 5 are illustrated in Figure 2. We can see that in V 3 the approximation is very poor since the frequencies are cut off excessively by the projection P 3 .If J is taken to be too large, the noise in the function g m is not damped enough by P J , and thus the high frequencies of g m are so extremely magnified that they destroy the approximated solution.The approximation parameter J 4 seems to be the optimal choice for this example.
In Figure 3 we display the exact solution, its approximation, and corresponding errors for k 5 and 100, respectively.We see that the proposed method is useful for different wave number k.
Figure 4 shows that the proposed method for the Cauchy problem for the modified Helmholtz equation is also effective.

Figure 1 :
Figure 1: a Exact solution u 1, • ; b unregularized solution reconstructed from g m for x 1.

Lemma 3 . 4 .
Let the function f λ : 0, a → R be given by c ∈ R and positive constants a < 1, b, and d.Then for the inverse function f −1 λ , one has

3 . 25 due to 3 . 21 .Theorem 3 . 5 .
Now, summarizing above inference process, we obtain the main result of the present paper.For s ≥ r, suppose that conditions 1.8 and 1.10 hold.If one takes J * * J * , 3.26

Figure 3 :
Figure 3: a and b correspond to k 5 and k 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 :
Figure 4: a Exact solution v 1, • ; b unregularized solution reconstructed from g m for x 1; c regularized solution reconstructed from g m for x 1 and J † 4; d the difference between the regularization and the exact solution.

Example 4 .
1. Take n 2 and g y e −y 2 ∈ S R 2 , where y y 1 , y 2 and S R 2 denotes the Schwartz function space.