Communications in Numerical Analysis

Volume 2013 (2013), Article ID cna-00178, 12 Pages

doi: 10.5899/2013/cna-00178

Research Article

Numerical solution of the helmholtz equation for the superellipsoid via the galerkin method

Yajni Warnapala1 *, Hy Dinh1

1Roger Williams University, Department of Mathematics, One Old Ferry Road, Bristol, RI 02809

* Corresponding author. Email address:; Tel. 401-254-3097.

Received: 19 November 2012; Accepted: 04 January 2013

Copyright © 2013 Yajni Warnapala and Hy Dinh. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The objective of this work was to find the numerical solution of the Dirichlet problem for the Helmholtz equation for a smooth superellipsoid. The superellipsoid is a shape that is controlled by two parameters. There are some numerical issues in this type of an analysis; any integration method is affected by the wave number k, because of the oscillatory behavior of the fundamental solution. In this case we could only obtain good numerical results for super ellipsoids that were more shaped like super cones, which is a narrow range of super ellipsoids. The formula for these shapes was: $x=cos(x)sin(y)^{n},y=sin(x)sin(y)^{n},z=cos(y)$ where $n$ varied from 0.5 to 4. The Helmholtz equation, which is the modified wave equation, is used in many scattering problems. This project was funded by NASA RI Space Grant for testing of the Dirichlet boundary condition for the shape of the superellipsoid. One practical value of all these computations can be getting a shape for the engine nacelles in a ray tracing the space shuttle. We are researching the feasibility of obtaining good convergence results for the superellipsoid surface. It was our view that smaller and lighter wave numbers would reduce computational costs associated with obtaining Galerkin coefficients. In addition, we hoped to significantly reduce the number of terms in the infinite series needed to modify the original integral equation, all of which were achieved in the analysis of the superellipsoid in a finite range. We used the Green's theorem to solve the integral equation for the boundary of the surface. Previously, multiple surfaces were used to test this method, such as the sphere, ellipsoid, and perturbation of the sphere, pseudosphere and the oval of Cassini Lin and Warnapala [9], Warnapala and Morgan [10].

Keywords: Helmholtz Equation; Galerkin Method; Superellipsoid Mathematics Subject Classifications (2000); 45B05; 65R10

1 Introduction

Figure 1: Superellipsoids for varying degrees of $n$

The main objective of this paper is to solve a boundary value problem for the Helmholtz equation. The Helmholtz equation is given by

\Delta u+k^{2}u=0,\ \text{Im }k\geq 0,

where $k$ is the wave number. Boundary value problems have being used for solving the Helmholtz equation, but this approach is less popular than the finite element method and the finite difference method. To overcome the non-uniqueness problem arising in integral equations for the exterior boundary-value problems for the Helmholtz's equation, Jones [5] suggested adding a series of outgoing waves to the free-space fundamental solution. Jost used this method for the Maxwell equations of electromagnetic scattering for the sphere with an explicit coefficient choice [6]. Here we use Jones modified integral equation approach, where we solve the exterior Dirichlet problem for the modified integral equation, using the same global Galerkin method used by Lin [8]. In this paper we looked at specifically the superellipsoid region, a versatile primitive which is controlled by two parameters. In this case we noted that there are numerical issues with very small and very large parameters. The shapes that we looked at were similar to the 3rd and the 4th shape from the left in (fig.1). Up to date there are no numerical results for the superellipsoid for the Helmholtz equation with the Dirichlet boundary condition.

2 Definitions

It can be shown that (Colton and Kress) that the Helmholtz integral equation is uniquely solvable if $k$ is not an eigenvalue for the corresponding interior Neumann problem. Therefore it is necessary to develop a method which is uniquely solvable for all frequencies $k.

Let $S$ be a closed bounded surface in $\Re ^{3}$ and assume it belongs to the class of $C^{2}$. Let $D_{-}$, $\ D_{+},$ denote the interior and exterior of $S$ respectively. The exterior Dirichlet problem for the Helmholtz's equation is given by

\begin{aligned} \Delta u(A)+k^{2}u(A) &=0, ~A=(x,y,z)\in D_{+},\text{Im }k\geq 0 \\ u(p) &=f(p),\ p\in S, \notag \end{aligned}


with $f$ a given function and $u$ satisfying the Sommerfeld radiation condition:

u=O(\frac{1}{r}),(\frac{\partial }{\partial r}-ik)u=o(\frac{1}{r})\text{ as } r=\left\vert A\right\vert \rightarrow \infty.


2.1 Theoretical Framework of the Boundary Value Problems

The exterior Dirichlet problem will be written as an integral equation. We represented the solution as a modified double layer potential, based on the modified fundamental solution. (See [4]).

u(A)=\int_{S}u(q)\frac{\partial (\frac{e^{ikr}}{4\pi r}+\chi (A,q))}{\partial \nu _{q}}d\sigma _{q}\text{with }A\in D_{+}\text{ whenever }=\left\vert A-q\right\vert.


The series of radiating waves is given by

\chi (A,q)=ik\sum_{n=0}^{\infty}\sum_{m=-n}^{n}a_{nm}h_{n}^{(1)}(k\left\vert A\right\vert )Y_{n}^{m}(\frac{A}{\left\vert A\right\vert })h_{n}^{(1)}(k\left\vert q\right\vert )\overline{Y_{n}^{m}}(\frac{q}{\left\vert q\right\vert }).


This addition of the infinite series to the fundamental solution is in order to remove the singularity that occurs when $A=q$.

Here $h_{n}^{(1)}$ denote the spherical Hankel function of the first kind and of order $n$, $Y_{n}^{m},$ $n=-m,...m$ are the linearly independent spherical harmonics of order $m$ given by $Y_{n}^{m}(\phi ,\theta )=(\dfrac{1}{2\pi}(m+\frac{1}{2})\dfrac{(m-n)!}{(m+n)!})^{\frac{1}{2}}p_{n}^{m}(\cos\theta )e^{im\phi}.$

As in [4] here we assume that $D_{-}$ (superellipsoid) to be a connected domain containing the origin and we choose a ball $B$ of radius $R$ and center at the origin such that $\overset{\_}{B}\subset D_{-}$. On the coefficients $a_{nm}$ we imposed the condition that the series $\chi (p,q)$ is uniformly convergent in $p$ and in $q$ in any region $\left\vert p\right\vert ,\left\vert q\right\vert \geq R+\epsilon $, $\epsilon >0$, and that the series can be two times differentiated term by term with respect to any of the variables with the resulting series being uniformly convergent. We also assumed that the series $\chi $ is a solution to the Helmholtz equation satisfying the Sommerfeld radiation condition for $\left\vert p\right\vert ,\left\vert q\right\vert >R.$ By letting $A$ tend to a point $ p\in S$, we obtain the following integral equation based on the Fredholm equations of the second kind

-2\pi \mu (p)+\int_{S}\mu (q)\frac{\partial \Psi (p,q)}{\partial \nu _{q}}d\sigma _{q}=-4\pi f(p),p\in S



\Psi =\frac{-e^{ikr_{qp}}}{r}-4\pi \chi (p,q).

We denote the above integral equations by

-2\pi \mu +K\mu =-4\pi f,


where in the Dirichlet case

K\mu (p)=\int_{S}\mu (q)\frac{\partial }{\partial \nu _{q}}(\frac{-e^{ikr_{qp}}}{r}-4\pi \chi (p,q))d\sigma _{q}.

By the assumptions on the series $\chi (p,q)$ the kernel $\frac{\partial \chi (p,q)}{\partial \nu _{q}}$ is continuous on $S\times S$, and hence $K$ is compact from $C(S)$ to $C(S)$ and $L^{2}(S)$ to $L^{2}(S).$

Kleinman and Roach [7] gave an explicit form of the coefficient $a_{nm}$ that minimizes the upper bound on the spectral radius (see [6]). If $B$ is the exterior of a sphere radius $R$ with center at the origin then the optimal coefficient for the Dirichlet problem was given by

a_{nm}=-\frac{1}{2}(\frac{j_{n}(kR)}{h_{n}^{(1)}(kR)}+\frac{j_{n}^{\prime}(kR)}{h_{n}^{(1){}^{\prime }}(kR)})\text{ for }n=0,1,2\text{ }...\text{ and }m=-n,...\text{ }n.

This choice of the coefficient minimizes the condition number, and (2.6) is uniquely solvable for the superellipsoid. This coefficient was given for spherical regions. Also the coefficient choice of $a_{nm}=-\dfrac{1}{2}(\dfrac{j_{n}(kR)}{h_{n}^{(1)}(kR)})$ was also considered but did not give good convergence results.

3 Smoothness of the integral operator $K$

Smoothness results of the double layer operator was proven by Lin [8]. We know that the series $\chi $ can be differentiated term by term with respect to any of the variables and that the resulting series is uniformly convergent. So the second derivative of the series is continuous on $\Re ^{3}\backslash B$ where $B=\left\{ x:\left\vert x\right\vert \leq R\right\} $. Furthermore the series $\chi $ is a solution to the Helmholtz equation satisfying the Sommerfeld radiation condition for $\left\vert x\right\vert ,\left\vert y\right\vert >R,$ when $B=\left\{ x:\left\vert x\right\vert \leq R\right\} $ is contained in $D.$

By (Theorem 3.5 [3]) any two times continuously differentiable solution of the Helmholtz's equation is analytic, and analytic functions are infinitely differentiable. So the series $\chi (p,q)$ is infinitely differentiable with respect to any of the variables $p,q$. Furthermore it is easy to see that if $\mu $ is bounded and integrable and $S\in C^{l}$ then $\int_{S}\chi (p,q)\mu (q)d\sigma _{q}\in C^{l}(S)$ and $\int_{S}\dfrac{\partial \chi (p,q)}{\partial \nu _{q}}\mu (q)d\sigma _{q}\in C^{l}(S).$

4 The Framework of the Galerkin Method

The variable of integration in (2.6) was changed, converting it to a new integral equation defined on the unit sphere. The Galerkin method was applied to this new equation,using spherical polynomials to define the approximating subspaces. $m:U\rightarrow_{onto}^{1-1}S,$ where $m$ is at least differentiable.

By changing the variable of integration on (2.6) we obtained the new equation over $U,$

-2\pi \widehat{\mu }+\widehat{K}\widehat{\mu }=-4\pi \widehat{f},\widehat{f}\in C(U).


The notation " $\symbol{94}$" denotes the change of variable from $S$ to $U.$ The operator $(-2\pi +\widehat{K})^{-1}$ exists and is bounded on $C(U)$ and $L^{2}(U).$ Let $X=L^{2}(U),\alpha =-2\pi ,$ and let an approximating subspace of spherical polynomials of degree $\leq N$ be denoted by $X_{N}$. The dimension of $X_{N}$ is $d_{N}=(N+1)^{2}:$ and we let $\left\{ h_{1},...h_{d}\right\}$ denote the basis of spherical harmonics Galerkin's method for solving (4.7) for the Dirichlet boundary conditions is given by

(-2\pi +P_{N}\widehat{K})\widehat{\mu _{N}}=-4\pi P_{N}\widehat{f.}


The solution is given by

\widehat{\mu _{N}}=\sum_{j=1}^{d}\alpha _{j}h_{j}

-2\pi \alpha _{i}(h_{i},h_{i})+\sum_{j=1}^{d}\alpha _{j}(\widehat{K} h_{j},h_{i})=-4\pi (\widehat{f},h_{i}),i=1,...\text{ }d.


The convergence of $\mu _{N}$ to $\mu $ in $L^{2}(S)$ is straightforward We know from previous literature that $P_{N}\widehat{\mu }\rightarrow \widehat{\mu }$ for all $\widehat{\mu }\in L^{2}(U)$. From standard results it follows that $\left\Vert \widehat{K}-P_{N}\widehat{K}\right\Vert \rightarrow 0$ and we can obtain the desired convergence.

(Also see [2]).Using the smoothness results of the integral operator $K$ from section III, and following the same proof as in [1], we can prove the following theorems.

Theorem 4.1 Assume that $f\in C^{l,\lambda }(S),$ $S\in C^{l+1,\lambda }$ $(S\in C^{2}$ for $l=0)$ and that the mapping $m$ satisfies (4.7) for some $l\geq 0.$ Then for all sufficiently large $N,$ the inverses $(-2\pi +P_{N}\widehat{K})^{-1}$ exist and are uniformly bounded and $\left\Vert \mu -\mu _{N}\right\Vert \leq \dfrac{c}{N^{l+\lambda ^{\prime }}}$ where $0< \lambda ^{{}^{\prime }}< \lambda $ is arbitrary. The constant $c$ depends on $l,\mu ,$ and $\lambda ^{\prime }.$

Convergence in $C(U).$ To prove uniform convergence of $\widehat{\mu _{N}}$ to $\widehat{\mu }$ is slightly more difficult. The main problem is that there are $\widehat{\mu }$ in $C(U)$ for which $P_{N}\widehat{\mu }$ does not converge to $\widehat{\mu .}$ Convergence for all $\widehat{\mu }$ would imply uniform boundedness of $\left\Vert P_{N}\right\Vert.$

Theorem 4.2 Assume that $S\in C^{2}$ and that $m$ satisfies (4.7) with $l=0.$ Then considering $\widehat{K}$ as an operator on $C(U),$

\left\Vert \widehat{K}-P_{N}\widehat{K}\right\Vert \rightarrow 0\text{ as }N\rightarrow \infty .


This implies the existence and uniform boundedness on $C(U)$ of $(-2\pi +P_{N}\widehat{K})^{-1}$ for all sufficiently large $N.$ Let $(-2\pi +\widehat{K})\widehat{\mu }=-4\pi \widehat{f}$ and $(-2\pi +P_{N}\widehat{K}) \widehat{\mu _{N}}=-4\pi P_{N}\widehat{f}.$ If $f\in C^{0,\lambda }(S),$ $\lambda >\frac{1}{2},$ then $\widehat{\mu _{N}}$ converges uniformly to $\widehat{\mu }.$ Moreover, if $S\in C^{l+1,\lambda }(S\in C^{2},$ for $l=0)$ and $f\in C^{l,\lambda }(S),$ $l+\lambda >\frac{1}{2},$ then $\left\Vert \mu -\mu _{N}\right\Vert _{\infty }\leq \dfrac{C}{N^{l+\lambda {}^{\prime }- \frac{1}{2}}}$ with $0< \lambda ^{\prime }< \lambda.$ The constant $c$ depends on $f,l,\lambda ^{\prime }.$

4.1 The approximation of true solutions for the Dirichlet problem

Given $\mu _{N}$ an approximate solution of (2.6), we defined the approximate solution $u_{N}$ of (2.1) using the integral (2.3).

u_{N}(A)=\int_{S}\mu _{N}(q)\frac{\partial }{\partial \upsilon _{q}}(\frac{e^{ikr_{qA}}}{4\pi r_{qA}}+\chi (A,q))d\sigma _{q},A\in D_{+}.


To show the convergence of $u_{N}(A),$ we used the following lemma.

Lemma 4.1

_{A\overset{\sup }{\in }K}\int_{S}\left\vert \frac{\partial }{\partial \nu_{q}}(\frac{e^{ikr_{qA}}}{4\pi r_{qA}}+\chi (A,q))\right\vert d\sigma_{q}< \infty ,


where $K$ is any compact subset of $D,$ from Warnapala and Morgan [10].

4.2 Implementation of the Galerkin Method for the Dirichlet Problem

The coefficients $(\widehat{K}h_{j},h_{i})$ are fourfold integrals with a singular integrand. Because the Galerkin coefficients $(\widehat{K} h_{j},h_{i})$ depends only on the surface $S,$ we calculated them separately for $N\leq N_{\max }.$ The following derivation was done for the Dirichlet problem. To decrease the effect of the singularity in computing $\widehat{K} h_{j}(\widehat{p})$ in the Dirichlet case, we used the identity

\int_{S}-\frac{\partial }{\partial \nu _{q}}\frac{1}{r_{qp}}d\sigma_{q}=2\pi ,p\in S\text{ where }r_{qp}=\left\vert p-q\right\vert ,

to write

\int_{S}-\frac{\partial }{\partial \nu _{q}}\frac{1}{r_{qp}}d\sigma_{q}=2\pi ,p\in S\text{ where }r_{qp}=\left\vert p-q\right\vert,\int_{U}-(h_{j}(\widehat{q})-h_{j}(\widehat{p})\widehat{\frac{\partial }{\partial \nu _{q}}\frac{1}{r_{qp}}}\left\vert J(\widehat{q})\right\vert d\sigma _{\widehat{q}}.

The integrands are bounded at $\widehat{q}=\widehat{p}$, where $J(\widehat{q})$ is the Jacobian.

5 Numerical Examples/Experimental Surfaces

In this section, several numerical examples are presented. The true solutions is given by


The parametric equation for the superellipsoid is given by $f(x,y)=(\cos (x)\sin (y)^{n},$ $\sin (x)\sin (y)^{n},cos(y))$ where $n$ varied from 0.5 to 4. The surface area of the superellipsoid has no closed form, thus one cannot give a specific formula but the volume is given by $\dfrac{8(\Gamma (1+\frac{1}{n}))^{3}}{\Gamma (1+\frac{1}{n})}$. For analysis it is important to realize that the superellipsoid is simply connected and is of infinite extent.

In our tables $INN$ are the interior nodes for calculating $\widehat{K} h_{j}, $ $EXN$ are the exterior nodes needed for calculating $(\widehat{K} h_{j,}h_{i})$ and $NINT$ are the nodes for calculating $u_{N}$. $N$ denotes the degree of the approximate spherical harmonics, recall that the number $d$ of basis functions equals to $(N+1)^{2}.$ In most cases we only added a few terms from the series. According to Jones [5] this is sufficient to remove the corresponding interior Dirichlet eigenvalues and obtain unique solutions at the same time.

Figure 2: The superellipsoid for $n=4$

Figure 3: The superellipsoid for $n=3$

Figure 4: The superellipsoid for $n=0.5$

The following Numerical results were obtained for the Dirichlet problem.

Table 1

$k=1,$ $n=3,$ $N=7,$ $INN=32,$ $EXN=16,$ true solution $u_{1}$


absolute error









In all the examples we used the coefficient $a_{nm}$.

Table 2

$k=1,$ $n=4$, $N=7,$ $INN=32,$ $EXN=16$, true solution $u_{1}$


absolute error









When we made the $k$ value smaller, we obtained the following Table 3(a).

Table 3(a)

$k=0.01,$ $n=3$, $N=7,$ $INN=32,$ $EXN=16,$

2nd boundary condition true solution $u_{1}$


absolute error









As you can see from the above Table 3(a) the results were mostly worse compared to the Table 1. Now we decreased the $n$ value and made the $k$ value an eigenvalue of the corresponding interior Neumann problem, and obtained Table 3(b).

Table 3(b)

$k=4.493409,$ $n=1.4$, $N=7,$ $INN=32,$ $EXN=16,$

number of terms $10$, true solution $u_{1}\medskip $


absolute error









From Table 3(b) it is evident that our method is successful as we get good convergence results for the eigenvalue of the interior Neumann problem. We also added more than five terms and still obtained good results. But as more terms and increasing of integration nodes, increases the CPU time considerably (this will be discussed later more extensively), we added only a few terms, only five in other cases of superellipsoid. In the next tables we changed the $n$ values.

Table 4

$k=0.001,$ $n=3$, $N=7,$ $INN=32,$ $EXN\ =16,$ true solution $u_{1}$


absolute error









From Table 1, 2, 3(a), 3(b) and Table 4, we see that for the points away from the boundary there is much greater accuracy than for points near the boundary. This is because the integrand is more singular at points near the boundary.


$k=0.01,$ $n=1.4$, $N=7,$ $INN=32,$ $EXN=16$, true solution $u_{1}$


absolute error









From Table 5 we see that to obtain similar accuracy as in the previous tables we might need to decrease the n value. This is due to the following fact: the superellipsoidal shape looks more spherical when the $n$ value decreases.

Remark 5.1 We picked $EXN< INN$, because the integrand of $(h_{i},\widehat{K}h_{j})$ is smoother than the integrand of $\widehat{K} h_{j}.$ We also picked $EXN\geq (N+1).$

More numerical data for varying values of $n.$

Table 6

$k=0.01,$ $n=2,$ $N=7,$ $INN=32,$ $EXN=16$, true solution $u_{1}$


absolute error









Table 7

$k=0.01,$ $n=1.8$, $N=7,$ $INN=32,$ $EXN=16,$ true solution $u_{1}$


absolute error









Table 8

$k=0.01,$ $n=1.5,$ $N=7,$ $INN=32,$ $EXN=16$, true solution $u_{1}$


absolute error









When $n$ is changed$,$ we obtained new shapes of superellipsoids. When we increased the nodes we obtained better results for relatively smaller $n$ values.

The (Fig.5) shows that further the points are from the boundary the better the convergence results are for this problem, which is the expected result.

Figure 5: Further the points are from the boundary better the convergence results

Remark 5.2 Few terms from the infinite series were added in all our numerical experiments. This is because in numerical calculations it is inefficient to add the full series. So we allow only a finite number of the coefficients $a_{nm}$ to be different from zero. According to Jones [5], this is sufficient to ensure uniqueness for the modified integral equations in a finite range of wave numbers $k.$ In practical applications, one is usually concerned with a finite range of $k$ so this is not a serious draw back. In order to use a large amount of nodes we need a considerably high amount of CPU time. From the above examples, we see that the error is effected by the boundary $S,$ $INN,$ $EXN$, boundary data and $k$ and $n$. As the value of the $n$ decreases the ellipsoidal nature of the superellipsoid for the Dirichlet problem increases and the rate of convergence decreases. If we want to obtain more accuracy, we must increase the number of integration nodes for calculating the Galerkin coefficients $(\widehat{K}h_{j},h_{i})$. The cost of calculating the Galerkin coefficients is high. Some of the increased cost comes from the complex number calculations, which is an intrinsic property of the Helmholtz equation. Furthermore any integration method is affected by $k$, due to the oscillatory behavior of the fundamental solution $\dfrac{e^{ikr}}{r}.$ Also the CPU time depends on the number of terms added from the series.

In order to eliminate more interior Neumann eigenvalues we need a more powerful computer which would decrease the CPU time considerably. We also see that for the shapes of superellipsoid with small $n$ values the convergence results are extremely good, even though the coefficient choice that was used was originally designed for spherical regions.


We would like to thank Dr. Rutherford and Mrs. Kennedy for their technical support in getting this paper ready. The author is grateful to the anonymous referees for their comments which substantially improved the quality of this paper.


  1. K. E. Atkinson, The numerical solution of the Laplace's equation in three dimensions, SIAM J. Numer. Anal, 19 (1982) 263-274.

  2. K. E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, (1998).

  3. David Colton, Rainer Kress, Integral equation methods in scattering theory, (1983).

  4. C. Criado, N. Alamo, Thomas rotation and Foucault pendulum under a simple unifying geometrical point of view, Int. J. Non-Linear Mech, 44 (2009) 923-927.

  5. D. S. Jones, Integral equations for the exterior acoustic problem, Q. J. Mech. Appl. Math, 27 (1974) 129-142.

  6. G. Jost, Integral equations with Modified Fundamental Solution in Time-Harmonic Electromagnetic Scattering, IMA J. Appl. Math, 40 (1988) 129-143.

  7. R. E. Kleinman, G. F. Roach, On modified Green functions in exterior problems for the Helmholtz equation, Proc. R. Soc. Lond. A, 383 (1982) 313-333.

  8. T. C. Lin, The numerical solution of Helmholtz equation using integral equations, Ph.D. thesis, University of Iowa, Iowa City, Iowa, (1982),

  9. T. C. Lin, Y. Warnapala-Yehiya, The Numerical Solution of the Exterior Dirichlet Problem for Helmholtz's Equation via Modified Green's Functions Approach, Comput. Math. Appl, 44 (2002) 1229-1248.

  10. Y. Warnapala, E. Morgan, The Numerical Solution of the Exterior Dirichlet Problem for Helmholtz's Equation via Modified Green's Functions Approach for the Oval of Cassini, Far East J. Appl. Math, 34 (2009) 1-20.