Solving the Cauchy problem for the Helmholtz equation using cubic smoothing splines

We consider the Cauchy problem for the Helmholtz equation defined in a rectangular domain. The Cauchy data are prescribed on a part of the boundary and the aim is to find the solution in the entire domain. The problem occurs in applications related to acoustics and is illposed in the sense of Hadamard. In our work we consider regularizing the problem by introducing a bounded approximation of the second derivative by using Cubic smoothing splines. We derive a bound for the approximate derivative and show how to obtain stability estimates for the method. Numerical tests show that the method works well and can produce accurate results. We also demonstrate that the method can be extended to more complicated domains.


Introduction
The Cauchy problem for the Helmholtz equation arises in many areas of applied science in relation to wave propagation, vibration and electromagnetic wave scattering. The commonly considered waves being water, sound and electromagnetic waves. It is often used in the determination of acoustic cavities [5] or the detection of the source of acoustical noise [6,23]. Other applications include the description of underwater waves [10], the determination of the radiation field surrounding a source of radiation [21], the localization of a tumor in a human brain [13], the identification and location of vibratory sources [15,16] and the detection of surface vibrations from interior acoustical pressure [7].
The Cauchy problem for the Helmholtz equation is ill-posed, in the sense of Hadamard, see [9], and therefore regularization is needed. In the literature, different regularization methods have been applied to the Cauchy problem for the Helmholtz equation. In [25] the problem was treated in the frequency domain, and an a priori bound u(·, a) 2 was used to stabilize the problem. A similar approach was used in [14], where the problem was stabilized using Meyer wavelets, which are compactly supported in the frequency domain, and also in [12]. In [24], the problem was regularized using a modified Tikhonov method. We also refer to [1][2][3][4]19,26], and the references therein, for more regularization techniques that have been successfully applied to the problem.
In this paper, we develop a novel regularization method for the Cauchy problem for the Helmholtz equation. In our previous work [18] we have shown that approximating a derivative by a bounded approximation based on Cubic smoothing splines has a stabilizing effect for the inverse heat conduction problem. Here we apply the same technique to the Helmholtz equation and are able to derive stability estimates. Also, we implement the method numerically and demonstrate that the method works well. Finally, we extend the method to the case of a modified Helmholtz equation, u+αu = 0, where α := α(x, y) > 0 is a function. The modified Helmholtz equation appears in situations where the problem is solved in more complicated domains .

The Helmholtz equation and ill-posedness
We consider the following Cauchy problem for the Helmholtz equation, in the rectangular domain = [0, 1] × [0, a]: where k 2 represents the wave number and g(x) and η(x) are the Cauchy data, available on part of the boundary of the domain.
The solution to (2.1), computed by the method of separation of variables, is where ω j = ( j 2 π 2 − k 2 ) 1/2 is a frequency, and are the Fourier-Sine coefficients. Note that the frequencies ω j may be either real or complex. The series (2.2) is convergent for the case of exact Cauchy data [g, η]. If we instead have noisy data [g δ , η δ ], such that g δ − g 2 and η δ − η 2 ≤ δ, where δ is a bound on the error, then there is no guarantee that the solution exists. This is because the magnification factors, e.g. | cosh(ω j y)| and | sinh(ω j y)/ω j |, for a fixed y > 0, grows exponentially as j → ∞. Thus a small error in the Fourier coefficients, i.e. g j and η j , for large values of j, are blown-up and may destroy the solution. This means that the problem is severely ill-posed, and stabilization is needed [9].
We remark that while the Cauchy problem (2.1) has a unique solution the same is not necessarily true for the boundary value problem obtained by replacing the Cauchy data, at y = 0, by Dirichlet data u(x, a) = f (x) and u(x, 0) = g(x). This can easily be seen by considering the function, which satisfies the Helmholtz equation in the domain , with the wave number k 2 = π 2 + (π/a) 2 and zero Dirichlet data. Similar solutions can be constructed for the case of zero Neumann, or Robin boundary, conditions. For our work we avoid issues related to non-uniqueness of the solution by only considering small values of the wave-number, i.e. k 2 ≤ π 2 (1 + 1/a 2 ), see [1].

Stabilization by discretizing the x-axis
In this section we stabilize the Cauchy problem for the Helmholtz equation by discretizing the x-variable. We replace the second derivative ∂ 2 x by a matrix D 2 in order to regularize the problem. For this purpose we rewrite (2.1) as an initial-value problem with initial-boundary data and u(0, y) = u(1, y) = 0, for 0 ≤ y ≤ a. (2.6) We discretize (2.4) on a uniform grid 0 = x 1 < · · · < x n = 1. For simplicity we introduce an operator that samples the functions on the grid such that, e.g g = (g(x)) = (g(x 1 ), . . . , g(x n )) T .
By introducing semi-discrete representations of the solution and its derivative, i.e.
U (y) = (u(·, y)) and U y (y) = (u y (·, y)) , we obtain the initial value problem, with initial data U (0) = g and U y (0) = η . The solution to (2.7) can formally be written as In [8] it was observed that since for any matrix the D 2 2 is bounded. Thus the system of ordinary differential equations (2.7) is well-posed and the solution depends in a stable way on the used data. By controlling the accuracy of the approximation D 2 the degree of stability can be adjusted to be suitable for a problem with a specific noise level.

Regularization by cubic smoothing splines
In this section we demonstrate that the Cauchy Problem for the Helmholtz equation can be regularized by using a cubic smoothing spline. We show how the approximation of ∂ 2 x ≈ D 2 is achieved using a cubic smoothing spline and also develop and prove Lemmas on the stability of the resulting regularized solution.

Differentiation by cubic smoothing splines
We adopt the basic notation about splines and introduce the second derivative matrix D 2 λ , show its properties and derive its bound in the Euclidean norm. We use a uniform grid, with grid parameter h. We work in an L 2 setting and introduce the Sobolev space together with the standard norms and semi-norms as defined in [20]. Let y = (y 1 , y 2 , . . . , y n ) T be a vector of data values on the grid. The cubic smoothing spline that approximates y on the grid is defined as follows: Definition 3.1 Let y ∈ R n be a vector. The cubic smoothing spline s λ [y] is obtained by solving min where λ > 0 is the regularization parameter and |u| 2 denotes the semi-norm.
It is known that the solution to (3.10) is a natural cubic spline and the corresponding operator, mapping y ∈ R n onto s λ [y] ∈ W 2 ([0, 1]), is linear, see [20,22]. Thus we can introduce a matrix approximating the second derivative as follows: The following properties of the matrix are useful: [18]. A direct computation shows that Since both β k (x) and β k (x) are even functions, with respect to x = x k , then we have . Hence rearranging the indicies shows that (D 2 λ ) i j = (D 2 λ ) ji , and the matrix is symmetric with real entries.
where c 3 is a constant.
Proof Since the second derivative of a cubic spline u is piecewise linear, with respect to the grid , we have h (u ) 2 2 ≤ c 2 1 |u | 2 0 , where c 1 is a constant. Then from the definition of D 2 λ we have where is the smoothing cubic spline corresponding to the vector y ∈ R n . Let y = y − u + u , then we can expand y 2 2 as A consequence of the fact that u(x) solves the least squares problem (3.10), is that h(y − u ) T u = λ|u| 2 2 , see [18,Lemma 3.7]. Inserting into (3.15), we obtain the estimate By using (3.16) in (3.14) we obtain the desired bound.

Stability analysis
In this section we establish stability results, for the initial value problem (2.7), when the approximation ∂ 2 x ≈ D 2 λ is used. The proofs are based on the idea that, since D 2 λ is a real symmetric matrix, there exists a unitary matrix X such that where the eigenvalues μ j are real. The stability results are established by a sequence of Lemmas: where λ is the regularization parameter, and let U 1 and U 2 be two different solutions to (2.7), corresponding to Cauchy data [(g 1 ) , η ] and The system (2.7) simplifies to Since λ is a diagonal matrix, the problem (3.19) can be solved for each eigenvalue μ j individually, with known boundary conditions V (μ j , 0) and V y (μ j , 0) = 0. We identify two distinct cases: where we write |k 2 +μ j | to emphasize the fact that we take the square root of a positive number. The constants C 1 and C 2 are computed by using the boundary conditions and we obtain the solution (3.22) Hence this case does not contribute to the instability of the problem.
Case 2 If μ j + k 2 < 0, we have the general solution as Again C 1 and C 2 are computed by making use of the boundary conditions and the solution is Now, since k 2 is positive and μ j +k 2 is assumed to be negative, we see that |μ j +k 2 | ≤ D 2 λ 2 − k 2 = σ 2 . Thus, due to cosh being a monotonically increasing function on the interval [0, ∞), we have (3.25) Here the magnification factor cosh(σ y) grows as μ j → −∞ and thus these components contribute to the instability of the problem. Finally, since X is orthogonal, we have V 2 = U 1 − U 2 2 and λ , with λ as the regularization parameter, and let U 1 and U 2 be two different solutions to (2.7), corresponding to Cauchy data [g , (η 1 The proof is similar to that of the previous Lemma. Introduce a new variable V (μ j , y) = X T (U 1 − U 2 ), and note that V (0) = X T (g − g) = 0 and V y (0) = X T (η 1 − η 2 ) . The problem is solved for one eigenvalue at a time, with known boundary conditions V y (μ j , 0) and V (μ j , 0) = 0. Again, we have two cases and if μ j + k 2 > 0 the solution can be written (3.29) By making use of the fact that | sin(α)/α| ≤ 1 we obtain |V (μ j , y)| ≤ y|V y (μ j , 0)|. (3.30) Thus the solution components satisfying μ j + k 2 > 0 are stable.

Proposition 3.7
Suppose that U λ is the regularized solution to (2.7), with exact data [g, η], and that U δ λ is the regularized solution with noisy data [g δ , η δ ], where λ is the regularization parameter. Then, if D 2 = D 2 λ and σ = D 2 λ 2 − k 2 > 0, Proof Let U 1 be the regularized solution to (2.7) with data [g, η δ ]. We have by the triangle inequality The result then follows by applying Lemma 3.5 and Lemma 3.6 to the respective terms.
The above proposition means that the approximation D 2 ≈ D 2 λ does stabalize the computations.

Simulated numerical examples
The numerical implementation of the splines based approach consists of implementing the procedure for computing the matrix-vector product D 2 λ y, for a given vector y ∈ R n . With the procedure in place, we are able to solve the initial value problem (2.7) using For all numerical experiments we used a uniform grid of size n = 500 in the xvariable. Note that we do not need to specify a grid in the y-direction since this is dealt with implicitly by the Runge-Kutta code. The specific tests are detailed below.
In order to investigate the ill-posedness of the inverse problem we added normally distributed noise of size δ = 10 −3 to the Cauchy data giving us noisy vectors [g δ , η δ ]. The inverse problem was solved for a wide range of regularization parameters λ, i.e. 10 −7 ≤ λ ≤ 10 −5 . For each value of the regularization parameter we computed a regularized solution f δ λ (x) := u δ λ (x, a) and also the error f δ λ − f 2 . In Fig. 2 we illustrate how the error depends on the parameter λ. We see that there is a trade-off between enforcing stability, i.e. a large regularization parameter λ, and solving the Helmholtz equation accurately, i.e. computing the second derivative accurately which means a small parameter λ. The conclusion is that by selecting an appropriate value for λ we are able to find a good regularized solution.
In actual applications the exact solution to the problem is unknown and we cannot compute the error f δ λ − f 2 . Thus another method for finding a good regularization parameter is needed. The L-curve, originally introduced in [17], is a plot of the norm of the approximate solution versus the residual in a log−log scale. Since our algorithm  (top-right). The L-curve suggests the parameter λ = 2.5 · 10 −7 and the minimum error is obtained for λ = 6.3 · 10 −7 . We also show the solutions obtained by using the λ value predicted by the L-curve (middle-left) and the λ value which gives the smallest error (middle-right). Finally, we show the solution obtained using λ = 2 · 10 −8 (bottom-left), which adds too little regularization, and λ = 1 · 10 −4 (bottom-right), which results in a too smooth solution. In all cases both the exact solution f (x) (black/solid-curves) and the inverse solution f δ λ (x) (blue/dashed-curves) are shown uses the semi-norm |u| 2 as the penalty term it is natural to measure the solution norm by the second derivative, ( f δ λ ) ) 2 . In our implementation the second derivative is approximated using the standard, centered, second order accurate, finite difference quotient. In order to define a residual we recall that the purpose of the residual is to measure the missmatch between an approximate solution and the available data. For our tests we define the residual as follows: Let v(x, t) satisfy the Helmholtz equation,  Fig. 3 The optimal regularization parameter λ as a function of the noise level δ present in the available Cauchy data (left graph). Also the error in the computed solution, when the optimal λ is used, as a function of the noise level (right graph). Here a = 0.2 and k 2 = 12 with Dirichlet data v(x, a) = f δ λ (x) and Neumann data v y (0, x) = η δ (x). The residual is then defined as v(x, 0) − g δ 2 . Note that this means that the residual is zero for the exact solution and data. However the residual does not necessarily tend to zero as λ → 0. In our work we compute the function v(x, t) by a standard finite difference scheme, see [1,2] for details. The computed L-curve is displayed in Fig. 2. We pick the regularization parameter λ that correspond to the point closest to the lower-left corner of the L-curve where both the solution norm and the residual is small. We see that the solution, obtained by using the parameter λ suggested by the L-curve, is nearly optimal.
Test 2 In order to further investigate the properties of the smoothing spline regularization method we carry out tests and find the optimal value of λ, and also the corresponding errors, as a function of the noise level. For this test the analytic solution from Test 1 was used with a = 0.2 and k 2 = 12. In the experiment we picked a large range of noise levels, i.e. 10 −4 ≤ δ ≤ 10 −1 , and for each level of noise we determined the optimal value of λ, i.e. the one that minimize the error f δ λ(δ) − f 2 . The results are shown in Fig. 3. We see that a smaller amount of noise δ means that we can use a smaller regularization parameter λ. This is consistent with the stability results derived in Sect. 3.2.

A modified Helmholtz equation
In this Section we look into solving the Cauchy problem for a modified Helmholz equation in the rectangle where α := α(x, y) is a positive, real valued function defined on the domain . The goal is to demonstrate that our Cubic spline regularization method can be used for solving this more general case. The modified Helmholtz equation is important in applications since it models wave propagation in a material with variable density. Also, consider a sitution where a function v satisfies the Helmholtz equation, with a constant wave number k 2 , in a more complicated domain , as illustrated in Fig. 4. For simplicity of notation we identify R 2 , with the complex plane C. If φ : → is the conformal mapping that transforms the domain into the rectangle then u(z) = v(φ −1 (z)), satisfies is understood as the derivative of an inverse function of a complex variable z. Thus, solving the modified Helmholtz equation in the rectangular domain is equivalent to solving the ordinary Helmholtz equation in the more complicated domain . Note that it is essential that the domains and are topologically equivalent in the sense that corners are mapped onto corners and smooth boundary segments are mapped onto the sides of the rectangle. Otherwise we would introduce singularities in the derivative boundary conditions, which is something the algorithm could not handle.
In order to demonstrate that the cubic spline regularization method can be applied to the problem (5.35) we select a function α(x, y) = 7 + sin(x/3) cos(y/2) + 3x y 2 ,  Fig. 5. Test 3 For this test the numerical test solution illustrated in Fig. 5 was used. Noisy Cauchy data was created by adding normally distributed random noise, with variance 10 −3 to the exact data vectors. We then solved the inverse problem for a wide range of regularization parameter values λ. The results are shown in Fig. 6. Again we see that there is an optimal regularization parameter, that represents the appropriate compromise between stability and accuracy. By using this parameter value we are able to reconstruct the solution accurately.  6 We display the error f δ λ − f 2 , as a function of λ (top-left). The error is minimized for λ = 2 · 10 −7 and the corresponding regularized solution f δ λ (x) is also displayed (top-right,blue curve) together with the exact solution f (x) (top-right,black-curve). The ill-posedness of the problem is illustrated by plotting the solution corresponding to λ = 5 · 10 −8 (bottom-left) which adds to little regularization. We also display the solution corresponding to λ = 8 · 10 −6 (bottom-right) which means too much added stability. In both cases the regularized solution f δ λ (blue/dashed curves) and the exact solution f (black/solid curves) are shown We also remark that experiments with different noise levels and different parameter values show that our method works well for the modified Helmholtz equation. The problem is severely ill-posed but still it is relatively easy to find a value for λ that works well.

Concluding remarks
In this paper we have developed a regularization method for solving the Cauchy problem for the Helmholtz equation based on Cubic smoothing splines. In our work we have developed a stability theory for the discrete problem that is solved numerically. One of the results is a matrix approximation of the second derivative for which we have obtained a norm-bound. The stability results are comparable to those obtained by different authors using similar methods.
Numerical experiments show that the method works well and can produce accurate results. For the experiments we have also introduced a residual that allows us to find appropriate regularization parameters by using the L-curve method. In our tests the parameters obtained from the L-curve are close to the optimal ones. We remark that the residual definition is clearly not optimal. Since we have to solve a well-posed boundary value problem to compute the residual we will have numerical errors. Thus, for a fixed step size for the discretization, the residual will not tend to zero as the noise level tends to zero. Also, since the Cauchy problem for the Helmholtz equation is severely ill-posed the results are quite sensitive with respect to the regularization parameter λ. More work is needed to obtain a good, and practical, parameter selection strategy.
The derived stability results, see Sect. 3.2, essentially show that as long as exp(ac 3 λ −1/4 )δ → 0, as δ → 0, then the propagated data error satisfies U λ(δ) − U δ λ(δ) 2 → 0. One potential parameter choice rule is which means exp(ac 3 λ −1/4 )δ = δ 1/2 . Since, by this rule, λ(δ) → 0 as δ → 0, we apply less and less regularization as the noise level is reduced, and thus it is resonable to expect that the appxoximation error U λ − U 2 also tends to zero. A formal proof of convergence for the above parameter choice rule is something we intend to do in future work. For our work we have primarily considered rectangular domains but in actual applications it is often the case that the Helmholtz equation is valid in a more complicated domain. Thus we extend our method to also treat the modified Helmholtz equation, i.e u + αu = 0, where α(x, y) > 0 is a function. This allows us to use a conformal mapping, or orthogonal mapping, to transform the domain under consideration to a rectangle. In this paper we only have initial numerical tests but this is something we intend to develop further in the future.
In the future we intend to apply our method to similar problems and make more systematic comparisons with the work of others, e.g. [3,11,14]. We also hope to apply the method to medical problems with measured data.