On the Computationally Efficient Numerical Solution to the Helmholtz Equation

Named after Hermann L. F. von Helmholtz (1821-1894), Helmholtz equation has obtained application in many fields: investigation of acaustic phenomena in aeronautics, electromagnetic application, migration in 3-D geophysical application, among many other areas. As shown in [2], Helmholtz equation is used in weather prediction at the Met Office in UK. Inefficiency, that is the bottleneck in Numerical Weather Prediction, arise partly from solving of the Helmholtz equation. This study investigates the computationally efficient iterative method for solving the Helmholtz equation. We begin by analysing the condition for stability of Jacobi Iterative method using Von Neumann method. Finally, we conclude that Bi-Conjugate Gradient Stabilised Method is the most computationally efficient method.


Introduction
In Numerical Weather Prediction, the discretisation of governing equations leads to the Helmholtz equation, which is solved for the increment in pressure field Π [3].The 2-dimensional Helmholtz equation is given by; where ∇ 2 is the Laplacian in rectangular coordinates, ω is the wave number and Φ is the value of the pressure field calculated from the previous time step.Solving of 3-dimensional version of the Helmholtz equation alone takes about a third of the total forecast time [2].As a result, optimal and efficient methods are being sought for solving this equation.In this study, a simple sample example of the Helmholtz problem is solved using Gauss-Seidel, Jacobi, Successive-Over-Relaxation, Conjugate Gradient, Bi-Conjugate Gradient, Bi-Conjugate Gradient Stabilised, Quasi-Minimal Residual and Generalised Minimal Residual Methods.Use is made of MATLAB for computation.Comparison of these methods is done in regards to computational efficiency.Before that, we analyse the stability of the 2-dimensional Helmholtz equation.

Discretisation of the Helmholtz equation
Let Φ(x, y) = 0.This is to ease the analysis of stability.Therefore, Central difference approximation of (1) yields the following difference equation: To simplify further, we let ∆x = ∆y = h.This substitution and multiplication by h 2 leads to: Collecting like terms together in (2.3) leads to: Equation (2.4) is the discretised form of the Helmholtz equation (2.1).

Stability analysis using Von Neumann method
Suppose we use Jacobi iteration for solving a system of equations generated by (2.4).The recurrence relation for the Helmholtz equation (2.1) is then given by: Let us define the error at the k th iteration as follows: Since the exact solution satisfies the relation (3.1), so is the error.Hence, The error term is represented in the form A pq is an arbitrary constant.Following from (3.4), we have the following expressions: In the selfsame way, Substituting equations (3.7) and (3.8) into equation (3.3) yields: ξ is the propagating factor.For stability; Equation (3.12) is the condition for stability of the Helmholtz equation (2.1).It shows that the wave number, ω, and the spatial step-size, h, determine the stability of the numerical scheme.To guarantee stability, the wave number and the spatial step-sizes should be chosen such that equation (3.12) is satisfied.

A sample problem
To analyse efficiency of various iterative methods, we make use of the following example.
Let us seek a numerical solution to the Helmholtz equation with boundary condition Let Φ(x, y) = 1 4 (x + y), ∆x = ∆y = 0.25 and ω = 1.Furthermore, let this problem be solved on a unit square, 0 ≤ x, y ≤ 1. Central difference discretisation of this problem, as illustrated in section 2 above, yields a system of equations in the form That is; In the following section, the system of equations (4.4) is solved using various iterative methods.A brief outline of the principle behind each method is given.

Iterative schemes
In this section, we outline and briefly describe various iterative schemes for solving the system of equations (4.4).More study on these methods can be obtained in [5] and [1] from which this summary is obtained.5.1.Jacobi Iterative (JI) Method.It is given by, where D, L and U are the diagonal, lower and upper matrices, respectively, obtained by decomposition of matrix A. k is the number of iteration.

5.2.
Gauss-Seidel Iterative (GSI) Method.This is given by, 2) The unknowns are similar to those in subsection (5.1) above.
5.3.Successive-Over Relaxation (SOR) Method.This is a modification of the Gauss-Seidel method to accelerate its convergence.It is expressed as follows: 3) α is called the relaxation factor.It can be shown that convergence is obtained within the range 0 < α < 2. To accelerate convergence of an already convergent Gauss-Seidel Method, over-relaxation method is used, that is, by choosing α in the interval 1 < α < 2. 5.4.Conjugate Gradient (CG) Method.This scheme is given by; where v (k) is a search direction expressed in vector form.
5.5.Bi-Conjugate Gradient Stabilized (BiCGSTAB) Method.This method avoids the often irregular convergence patterns of the Conjugate Gradient Squared method.
where Q i is an ith degree polynomial describing a deepest descent.
5.6.Bi-Conjugate Gradient (BICG) Method.Unlike in the Conjugate Gradient method, the residuals are replaced with relations that are similar but based on A T instead of A .Two sequences of residuals are updated.
This ensures the bi-orthogonality relations: It is meant to overcome the irregularity of the convergence of BICG method.It solves the reduced tri-diagonal system in a least squares sense.It uses look-ahead techniques to avoid breakdowns that make it more robust than BICG.

Results
The solution of the sample 2-dimensional Helmholtz equation has been summarised in Table 1 and Figure 1 for the different methods used.Solving the Helmholtz equation iteratively, using different methods, shows that the iterations converge faster by Bi-Conjugate Gradient Stabilized Method followed by the Generalised Minimal Residual Method.The least convergent method is the Jacobi Iterative Method.While the Quasi-Minimal Residual Method converges relatively faster than the Jacobi Iterative Method, it takes longer CPU-time than all other methods.It is therefore not computationally efficient.The most computationally efficient method for solving the Helmholtz equation is therefore Bi-Conjugate Gradient Stabilized Method followed by Generalised Minimal Residual Method.

Figure 1 .
Figure 1.A graph showing the number of iterations and CPUtime of different methods

Figure 2 .
Figure 2. Visualization of the Helmholtz equation in MATLAB

Table 1 .
A summary of the number of iterations and CPU-Time for methods used 11The zero value in Table1implies a negligible value