Fourier-Hankel solution of the Robin problem for the Helmholtz equation in supershaped annular domains

Abstract

The Robin problem for the Helmholtz equation in normal-polar annuli is addressed by using a suitable Fourier-Hankel series technique. Attention is in particular focused on the wide class of domains whose boundaries are defined by the so-called superformula introduced by Gielis. A dedicated numerical procedure based on the computer algebra system Mathematica^© is developed in order to validate the proposed methodology. In this way, highly accurate approximations of the solution, featuring properties similar to the classical ones, are obtained.

Introduction

Many problems of mathematical physics and electromagnetics are related to the Laplacian:

• The wave equation $v_{tt}=a^2\mathrm{\Delta }v$;

• The heat propagation $v_t=\kappa \mathrm{\Delta }v$;

• The Laplace equation $\mathrm{\Delta }v=0$;

• The Poisson equation $\mathrm{\Delta }v=f$;

• The Helmholtz equation $\mathrm{\Delta }v+k^{2}v=0$;

• The Schrödinger equation $-\frac{h^2}{2m}\mathrm{\Delta }\psi +V\psi =E\psi$.

In recent papers [1–8], the classical Fourier projection method [9, 10] for solving boundary value problems (BVP s) for the Laplace and Helmholtz equations in canonical domains has been extended in order to address similar differential problems in simply connected starlike domains, whose boundaries may be regarded as an anisotropically stretched unit circle centered at the origin.

In this contribution, a suitable technique useful to compute the coefficients of the Fourier-Hankel expansion representing the solution of the Robin boundary value problem for the Helmholtz equation in complex annular domains is presented. In particular, the boundaries of the considered domains are supposed to be defined by the so-called Gielis formula [11]. Regular functions are assumed to describe the boundary values, but the proposed approach can be easily generalized in case of weakened hypotheses.

In order to verify and validate the developed methodology, a suitable numerical procedure based on the computer algebra system Mathematica^© has been adopted. By using such a procedure, a point-wise convergence of the Fourier-Hankel series representation of the solution has been observed in the regular points of the boundaries, with Gibbs-like phenomena potentially occurring in the quasi-cusped points. The obtained numerical results are in good agreement with theoretical findings by Carleson [12].

The Laplacian in stretched polar coordinates

Let us introduce in the real plane the usual polar coordinate system

$$\{\begin{array}{c}x=rcos\vartheta ,\hfill \\ y=rsin\vartheta ,\hfill \end{array}$$

(1)

and the polar equations

$$r=R_{\pm }(\vartheta ),$$

(2)

relevant to the boundaries of the supershaped annulus which is described by the following chain of inequalities:

$$R_{-}(\vartheta )\le r\le R_{+}(\vartheta ),$$

(3)

with $0\le \vartheta \le 2\pi$. In (2), $R_{\pm }(\vartheta )$ are assumed to be piece-wise $C^2$ functions satisfying the condition

$$R_{+}(\vartheta )>R_{-}(\vartheta )>0,0\le \vartheta \le 2\pi .$$

(4)

In this way, upon introducing the stretched radius ϱ such that

$$r=\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a},$$

(5)

with $b>a>0$, the considered annular domain can be readily obtained by taking $a\le \varrho \le b$.

Remark Note that in the stretched coordinate system ϱ, ϑ, the original domain is transformed into the circular annulus of radii a and b, respectively. Hence, in this system one can use classical techniques, including the eigenfunction method, to solve the Helmholtz equation [10].

Let us consider a piece-wise $C^2(\mathring{\mathcal{A}})$ function $v(x,y)=v(rcos\vartheta ,rsin\vartheta )=u(r,\vartheta )$ and the Laplace operator in polar coordinates

$$\mathrm{\Delta }u=\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\vartheta }^2}.$$

(6)

In the considered stretched coordinate system, Δ can be represented by setting

$$U(\varrho ,\vartheta )=u(\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a},\vartheta ).$$

(7)

In this way, by denoting $R_{\pm }(\vartheta )$ as $R_{\pm }$ for the sake of shortness, one can readily find:

$$\frac{\partial u}{\partial r}=\frac{b-a}{R_{+}-R_{-}}\frac{\partial U}{\partial \varrho },$$

(8)

$$\frac{{\partial }^{2}u}{\partial r^2}={\left(\frac{b-a}{R_{+}-R_{-}}\right)}^2\frac{{\partial }^{2}U}{\partial {\varrho }^2},$$

(9)

$$\begin{array}{c}\frac{{\partial }^{2}u}{\partial {\vartheta }^2}=\{\frac{(b-\rho )[2{\dot{R}}_{-}({\dot{R}}_{+}-{\dot{R}}_{-})-{\ddot{R}}_{-}(R_{+}-R_{-})]}{{(R_{+}-R_{-})}^2}\hfill \\ \phantom{\frac{{\partial }^{2}u}{\partial {\vartheta }^2}=}-\frac{(a-\rho )[2{\dot{R}}_{+}({\dot{R}}_{+}-{\dot{R}}_{-})-{\ddot{R}}_{+}(R_{+}-R_{-})]}{{(R_{+}-R_{-})}^2}\}\frac{\partial U}{\partial \varrho }\hfill \\ \phantom{\frac{{\partial }^{2}u}{\partial {\vartheta }^2}=}+{\left[\frac{(b-\varrho ){\dot{R}}_{-}-(a-\varrho ){\dot{R}}_{+}}{R_{+}-R_{-}}\right]}^2\frac{{\partial }^{2}U}{\partial {\varrho }^2}-2\frac{(b-\varrho ){\dot{R}}_{-}-(a-\varrho ){\dot{R}}_{+}}{R_{+}-R_{-}}\frac{{\partial }^{2}U}{\partial \varrho \partial \vartheta }+\frac{{\partial }^{2}U}{\partial {\vartheta }^2},\hfill \end{array}$$

(10)

where the dot superscript denotes the differentiation with respect to the angle ϑ. Substituting equations (8)-(10) into equation (6) finally yields

$$\begin{array}{rl}\mathrm{\Delta }u=& {\left(\frac{b-a}{R_{+}-R_{-}}\right)}^2(\{1+{\left[\frac{(b-\varrho ){\dot{R}}_{-}-(a-\varrho ){\dot{R}}_{+}}{(b-\varrho )R_{-}-(a-\varrho )R_{+}}\right]}^2\}\frac{{\partial }^{2}U}{\partial {\varrho }^2}\\ +\{\frac{(b-\rho )[(R_{-}-{\ddot{R}}_{-})(R_{+}-R_{-})+2{\dot{R}}_{-}({\dot{R}}_{+}-{\dot{R}}_{-})]}{{[(b-\varrho )R_{-}-(a-\varrho )R_{+}]}^2}\\ -\frac{(a-\rho )[(R_{+}-{\ddot{R}}_{+})(R_{+}-R_{-})+2{\dot{R}}_{+}({\dot{R}}_{+}-{\dot{R}}_{-})]}{{[(b-\varrho )R_{-}-(a-\varrho )R_{+}]}^2}\}\frac{\partial U}{\partial \varrho }\\ -2(R_{+}-R_{-})\frac{(b-\varrho ){\dot{R}}_{-}-(a-\varrho ){\dot{R}}_{+}}{{[(b-\varrho )R_{-}-(a-\varrho )R_{+}]}^2}\frac{{\partial }^{2}U}{\partial \varrho \partial \vartheta }\\ +{\left[\frac{R_{+}-R_{-}}{(b-\varrho )R_{-}-(a-\varrho )R_{+}}\right]}^2\frac{{\partial }^{2}U}{\partial {\vartheta }^2}).\end{array}$$

(11)

As it can be easily noticed, upon setting $R_{-}(\vartheta )=a=0$ and $R_{+}(\vartheta )=b=1$, the classical expression of the Laplacian in polar coordinates is recovered.

The Robin problem for the Helmholtz equation

Let us consider the interior Robin problem for the Helmholtz equation in a starlike annulus , whose boundaries ${\partial }_{\pm }\mathcal{A}$ are described by the polar equations $r=R_{\pm }(\vartheta )$ respectively:

$$\{\begin{array}{cc}\mathrm{\Delta }v(x,y)+k^{2}v(x,y)=0,\hfill & (x,y)\in \mathring{\mathcal{A}},\hfill \\ {\lambda }_{\pm }v(x,y)+{\gamma }_{\pm }\frac{\partial v}{\partial \nu }(x,y)=f_{\pm }(x,y),\hfill & (x,y)\in {\partial }_{\pm }\mathcal{A},\hfill \end{array}$$

(12)

where $k>0$ denotes the propagation constant, ${\hat{\nu }}_{\pm }={\hat{\nu }}_{\pm }(\vartheta )$ are the outward-pointing normal unit vectors to the domain boundaries ${\partial }_{\pm }\mathcal{A}$, respectively, and ${\lambda }_{\pm }$, ${\gamma }_{\pm }$ are given regular weighting coefficients.

Under the mentioned assumptions, one can prove the following theorem.

Theorem Let

$${\psi }_{\pm }(\vartheta )=\frac{d}{d\vartheta }ln{R}_{\pm }(\vartheta )=\frac{{\dot{R}}_{\pm }(\vartheta )}{R_{\pm }(\vartheta )},$$

(13)

and

$$\begin{array}{rl}{f}_{\pm }(R_{\pm }(\vartheta )cos\vartheta ,R_{\pm }(\vartheta )sin\vartheta )& =F_{\pm }(\vartheta )\\ =\frac{1}{\sqrt{1+{\psi }_{\pm }{(\vartheta )}^2}}\sum _{m=0}^{+\mathrm{\infty }}({\alpha }_m^{(\pm )}cosm\vartheta +{\beta }_m^{(\pm )}sinm\vartheta ),\end{array}$$

(14)

where

$$\left\{\begin{array}{c}{\alpha }_m^{(\pm )}\\ {\beta }_m^{(\pm )}\end{array}\right\}=\frac{{ϵ}_m}{2\pi }{\int }_0^{2\pi }{F}_{\pm }(\vartheta )\sqrt{1+{\psi }_{\pm }{(\vartheta )}^2}\left\{\begin{array}{c}cosm\vartheta \\ sinm\vartheta \end{array}\right\}d\vartheta ,$$

(15)

${ϵ}_m$ being the usual Neumann symbol. Then boundary value problem (12) for the Helmholtz equation admits a classical solution $v(x,y)\in L^2(\mathcal{A})$ such that the following Fourier-Hankel series expansion holds true:

$$\begin{array}{c}v(\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}cos\vartheta ,\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}sin\vartheta )\hfill \\ =U(\varrho ,\vartheta )\hfill \\ =\sum _{m=0}^{+\mathrm{\infty }}[H_m^{(1)}\left(k\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}\right)(A_{1,m}cosm\vartheta +B_{1,m}sinm\vartheta )\hfill \\ +H_m^{(2)}\left(k\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}\right)(A_{2,m}cosm\vartheta +B_{2,m}sinm\vartheta )].\hfill \end{array}$$

(16)

For each index m, define

$$\begin{array}{rl}\left[\begin{array}{c}{\xi }_{p,m}^{(\pm )}(\vartheta )\\ {\eta }_{p,m}^{(\pm )}(\vartheta )\end{array}\right]=& \left[\begin{array}{cc}cosm\vartheta & -sinm\vartheta \\ sinm\vartheta & cosm\vartheta \end{array}\right]\\ \cdot \left[\begin{array}{c}{\lambda }_{\pm }{H}_m^{(p)}(k{R}_{\pm }(\vartheta ))\sqrt{1+{\psi }_{\pm }{(\vartheta )}^2}\pm k{\gamma }_{\pm }{\dot{H}}_m^{(p)}(k{R}_{\pm }(\vartheta ))\\ \mp {\gamma }_{\pm }\frac{m}{R_{\pm }(\vartheta )}{H}_m^{(p)}(k{R}_{\pm }(\vartheta )){\psi }_{\pm }(\vartheta )\end{array}\right],\end{array}$$

(17)

with $H_m^{(p)}(\cdot )$ denoting the Hankel function of kind $p=1,2$ and order m. Hence, the coefficients $A_{p,m}$, $B_{p,m}$ in (16) can be determined by solving the infinite linear system

$$\sum _{m=0}^{+\mathrm{\infty }}\left[\begin{array}{cccc}{\mathrm{X}}_{1,1,n,m}^{(-)}& {\mathrm{Y}}_{1,1,n,m}^{(-)}& {\mathrm{X}}_{1,2,n,m}^{(-)}& {\mathrm{Y}}_{1,2,n,m}^{(-)}\\ {\mathrm{X}}_{2,1,n,m}^{(-)}& {\mathrm{Y}}_{2,1,n,m}^{(-)}& {\mathrm{X}}_{2,2,n,m}^{(-)}& {\mathrm{Y}}_{2,2,n,m}^{(-)}\\ {\mathrm{X}}_{1,1,n,m}^{(+)}& {\mathrm{Y}}_{1,1,n,m}^{(+)}& {\mathrm{X}}_{1,2,n,m}^{(+)}& {\mathrm{Y}}_{1,2,n,m}^{(+)}\\ {\mathrm{X}}_{2,1,n,m}^{(+)}& {\mathrm{Y}}_{2,1,n,m}^{(+)}& {\mathrm{X}}_{2,2,n,m}^{(+)}& {\mathrm{Y}}_{2,2,n,m}^{(+)}\end{array}\right]\cdot \left[\begin{array}{c}{A}_{1,m}\\ B_{1,m}\\ A_{2,m}\\ B_{2,m}\end{array}\right]=\left[\begin{array}{c}{\alpha }_n^{(-)}\\ {\beta }_n^{(-)}\\ {\alpha }_n^{(+)}\\ {\beta }_n^{(+)}\end{array}\right],$$

(18)

where

$${\mathrm{X}}_{\left\{\begin{array}{c}1\\ 2\end{array}\right\},p,n,m}^{(\pm )}=\frac{{ϵ}_n}{2\pi }{\int }_0^{2\pi }{\xi }_{p,m}^{(\pm )}(\vartheta )\left\{\begin{array}{c}cosn\vartheta \\ sinn\vartheta \end{array}\right\}d\vartheta ,$$

(19)

$${\mathrm{Y}}_{\left\{\begin{array}{c}1\\ 2\end{array}\right\},p,n,m}^{(\pm )}=\frac{{ϵ}_n}{2\pi }{\int }_0^{2\pi }{\eta }_{p,m}^{(\pm )}(\vartheta )\left\{\begin{array}{c}cosn\vartheta \\ sinn\vartheta \end{array}\right\}d\vartheta ,$$

(20)

with $p=1,2$ and $m,n\in {\mathbb{N}}_0:=\mathbb{N}\cup \{0\}$.

Proof Upon noting that in the stretched coordinate system ϱ, ϑ introduced in the x, y plane, the considered domain turns into the circular annulus of radii a and b, one can readily adopt the usual eigenfunction method [10] in combination with the separation of variables (with respect to r and ϑ). As a consequence, elementary solutions of the problem can be searched in the form

$$u(r,\vartheta )=U(\frac{b[r-R_{-}(\vartheta )]-a[r-R_{+}(\vartheta )]}{R_{+}(\vartheta )-R_{-}(\vartheta )},\vartheta )=\mathrm{P}(r)\mathrm{\Theta }(\vartheta ).$$

(21)

Substituting into the Helmholtz equation, one easily finds that the functions $\mathrm{P}(\cdot )$, $\mathrm{\Theta }(\cdot )$ must satisfy the ordinary differential equations

$$\frac{d^2\mathrm{\Theta }(\vartheta )}{d{\vartheta }^2}+{\mu }^2\mathrm{\Theta }(\vartheta )=0,$$

(22)

$$r^2\frac{d^2\mathrm{P}(r)}{d{r}^2}+r\frac{d\mathrm{P}(r)}{dr}+(k^2{r}^2-{\mu }^2)\mathrm{P}(r)=0,$$

(23)

respectively. The parameter μ is a separation constant whose choice is governed by the physical requirement that at any fixed point in the real plane the scalar field $u(r,\vartheta )$ must be single-valued. So, by setting $\mu =m\in {\mathbb{N}}_0$, one can easily find

$$\mathrm{\Theta }(\vartheta )=a_{m}cosm\vartheta +b_{m}sinm\vartheta ,$$

(24)

where $a_m,b_m\in \mathbb{C}$ denote arbitrary constants. The radial function $\mathrm{P}(\cdot )$ satisfying (23) can be readily expressed as follows:

$$\mathrm{P}(r)=c_m{H}_m^{(1)}(kr)+d_m{H}_m^{(2)}(kr),$$

(25)

with $c_m,d_m\in \mathbb{C}$. Therefore, the general solution of Robin problem (12) can be searched in the form

$$\begin{array}{rl}u(r,\vartheta )=& \sum _{m=0}^{+\mathrm{\infty }}[H_m^{(1)}(kr)(A_{1,m}cosm\vartheta +B_{1,m}sinm\vartheta )\\ +H_m^{(2)}(kr)(A_{2,m}cosm\vartheta +B_{2,m}sinm\vartheta )].\end{array}$$

(26)

Enforcing the Robin boundary condition yields

$$\begin{array}{rl}{F}_{\pm }(\vartheta )& ={\lambda }_{\pm }u(R_{\pm }(\vartheta ),\vartheta )+{\gamma }_{\pm }\frac{\partial u}{\partial \nu }(R_{\pm }(\vartheta ),\vartheta )\\ ={\lambda }_{\pm }u(R_{\pm }(\vartheta ),\vartheta )+{\gamma }_{\pm }\mathrm{\nabla }u(R_{\pm }(\vartheta ),\vartheta )\cdot {\hat{\nu }}_{\pm }(\vartheta ),\end{array}$$

(27)

where

$$\mathrm{\nabla }u(r,\vartheta )=\hat{r}\frac{\partial u(r,\vartheta )}{\partial r}+\hat{\vartheta }\frac{1}{r}\frac{\partial u(r,\vartheta )}{\partial \vartheta },$$

(28)

and

$${\hat{\nu }}_{\pm }(\vartheta )=\pm \frac{\hat{r}-{\psi }_{\pm }(\vartheta )\hat{\vartheta }}{\sqrt{1+{\psi }_{\pm }{(\vartheta )}^2}}.$$

(29)

Hence, combining equations above and using the classical Fourier projection method, equations (17)-(20) follow after some algebraic manipulations. □

It is worth noting that the derived expressions still hold under the assumption that $R_{\pm }(\vartheta )$ are piecewise continuous functions, and the boundary values are described by square integrable, not necessarily continuous, functions, so that the relevant Fourier coefficients ${\alpha }_m^{(\pm )}$, ${\beta }_m^{(\pm )}$ in equation (14) are finite quantities.

Numerical procedure

In the following numerical examples, let us assume, for the boundaries ${\partial }_{\pm }\mathcal{A}$ of the considered annulus, general polar equations of the type

$$R_{\pm }(\vartheta )={\left(\right|\frac{1}{d_x^{\pm }}cos\frac{k_x^{\pm }\vartheta }{4}{|}^{{\nu }_x^{\pm }}+|\frac{1}{d_y^{\pm }}sin\frac{k_y^{\pm }\vartheta }{4}{|}^{{\nu }_y^{\pm }})}^{-1/{\nu }_0^{\pm }},$$

(30)

as introduced by Gielis in [11]. Very different characteristic geometries, including ellipses, Lamé curves, ovals, and m-fold symmetric figures, are obtained by assuming suitable values of the parameters $k_x^{\pm }$, $k_y^{\pm }$, $d_x^{\pm }$, $d_y^{\pm }$, ${\nu }_x^{\pm }$, ${\nu }_y^{\pm }$, ${\nu }_0^{\pm }$ in (30). It is emphasized that almost all two-dimensional normal-polar annular domains can be described, or closely approximated, by the above-mentioned Gielis formula.

In order to assess the performance of the proposed methodology in terms of numerical accuracy and convergence rate, the relative boundary error has been evaluated as follows:

$$\begin{array}{rl}{e}_N=& \frac{\parallel {\lambda }_{-}{U}_N(a,\vartheta )+{\gamma }_{-}\frac{\partial U_N}{\partial \nu }(a,\vartheta )-F_{-}(\vartheta )\parallel }{\parallel F_{-}(\vartheta )\parallel }\\ +\frac{\parallel {\lambda }_{+}{U}_N(b,\vartheta )+{\gamma }_{+}\frac{\partial U_N}{\partial \nu }(b,\vartheta )-F_{+}(\vartheta )\parallel }{\parallel F_{+}(\vartheta )\parallel },\end{array}$$

(31)

with $\parallel \cdot \parallel$ being the usual $L^2$ norm, and where $U_N(\varrho ,\vartheta )$ denotes the partial sum of order N relevant to the Fourier-Hankel series expansion representing the solution of the boundary value problem for the Helmholtz equation, namely

$$\begin{array}{r}{U}_N(\varrho ,\vartheta )\\ =\sum _{m=0}^N[H_m^{(1)}\left(k\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}\right)(A_{1,m}cosm\vartheta +B_{1,m}sinm\vartheta )\\ +H_m^{(2)}\left(k\frac{(b-\varrho )R_{-}(\vartheta )-(a-\varrho )R_{+}(\vartheta )}{b-a}\right)(A_{2,m}cosm\vartheta +B_{2,m}sinm\vartheta )].\end{array}$$

(32)

Remark It is to be noticed that where the boundary values exhibit a rapidly oscillating behavior, the number N of terms in expansion (32) approximating the solution of the problem should be increased accordingly in order to achieve the desired numerical accuracy.

First example

By assuming in (30) $k_x^{\pm }=k_y^{\pm }=4$, $d_x^{-}=d_y^{-}=3/4$, $d_x^{+}=d_y^{+}=5/2$, ${\nu }_x^{\pm }={\nu }_y^{\pm }=6$, ${\nu }_0^{\pm }=10$, the annulus features a strip-like shape. Let $f_{-}(x,y)=-x^2+y+2i{y}^2$ and $f_{+}(x,y)=cos(x/2)+isin(x/2)$ be the functions describing the boundary values. Provided that the propagation constant is $k=1$, and ${\lambda }_{-}=1$, ${\lambda }_{+}=-1$, ${\gamma }_{-}=1/10$, ${\gamma }_{+}=-2$ are the weighting coefficients in the Robin condition, the relative boundary error $e_N$ as a function of the number N of terms in truncated series expansion (32) exhibits the behavior shown in Figure 1. As it appears from Figure 2, the selection of the expansion order $N=10$ leads to a very accurate Fourier-Hankel representation $v_N(x,y)$ of the solution (featuring $e_N<1\mathrm{\%}$). The spatial distribution of $v_N(x,y)$ is shown in Figure 3, whereas the logarithmic magnitude of the relevant expansion coefficients $A_{p,m}$ and $B_{p,m}$ ($p=1,2$) is plotted in Figure 4.

Figure 1

Relative boundary error ${\mathit{e}}_{\mathit{N}}$ as a function of order N of the truncated Fourier-Hankel series expansion representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{5}\mathbf{/}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{6}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{\pm }}\mathbf{=}\mathbf{10}$ .

Figure 2

Boundary behavior along ${\mathit{\partial }}_{\mathbf{-}}\mathcal{A}$ (a) and ${\mathit{\partial }}_{\mathbf{+}}\mathcal{A}$ (b) of the partial sum ${\mathit{U}}_{\mathit{N}}\mathbf{(}\mathit{\varrho }\mathbf{,}\mathit{\vartheta }\mathbf{)}$ of order $\mathit{N}\mathbf{=}\mathbf{10}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{5}\mathbf{/}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{6}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{\pm }}\mathbf{=}\mathbf{10}$ .

Figure 3

Spatial distribution of the Fourier-Hankel series expansion ${\mathit{v}}_{\mathit{N}}\mathbf{(}\mathit{x}\mathbf{,}\mathit{y}\mathbf{)}$ of order $\mathit{N}\mathbf{=}\mathbf{10}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{5}\mathbf{/}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{6}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{\pm }}\mathbf{=}\mathbf{10}$ .

Figure 4

Logarithmic magnitude of the coefficients ${\mathit{A}}_{\mathit{p}\mathbf{,}\mathit{m}}$ and ${\mathit{B}}_{\mathit{p}\mathbf{,}\mathit{m}}$ ( $\mathit{p}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{2}$ ) relevant to the Fourier-Hankel series expansion ${\mathit{v}}_{\mathit{N}}\mathbf{(}\mathit{x}\mathbf{,}\mathit{y}\mathbf{)}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{5}\mathbf{/}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{\pm }}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{\pm }}\mathbf{=}\mathbf{6}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{\pm }}\mathbf{=}\mathbf{10}$ .

Second example

In the second numerical example, we turn to the consideration of the class of annuli having one or both boundaries featuring a polygonal contour. In this respect, it is not difficult to show that the general k-sided convex regular polygon can be readily described by the following specialized version of Gielis’ formula [13]:

$$R_k(\vartheta )=\underset{\nu \to +\mathrm{\infty }}{lim}{\left[\right|\frac{1}{d}cos\frac{k\vartheta }{4}{|}^{2(1-\nu {log}_{2}cos\frac{\pi }{k})}+|\frac{1}{d}sin\frac{k\vartheta }{4}{|}^{2(1-\nu {log}_{2}cos\frac{\pi }{k})}]}^{-1/\nu }.$$

(33)

In this way, the methodology detailed in the previous section can be used straightforwardly. In particular, upon assuming in (30) $k_x^{-}=k_y^{-}=2$, $d_x^{-}=1/2$, $d_y^{-}=3/4$, ${\nu }_x^{-}=2$, ${\nu }_y^{-}={\nu }_0^{-}=3$, as well as $k_x^{+}=k_y^{+}=k^{+}=5$, $d_x^{+}=d_y^{+}=9/4$, and ${\nu }_x^{+}={\nu }_y^{+}=2(1-{\nu }_0^{+}{log}_{2}cos\frac{\pi }{k^{+}})$, with ${\nu }_0^{+}\to +\mathrm{\infty }$, the annulus may be regarded as the result of the Boolean subtraction of an ovaloid from a regular pentagon. Let $f_{-}(x,y)=1/(i{e}^{x^2+y}+y^2)$ and $f_{+}(x,y)=1$ be the functions describing the boundary values along ${\partial }_{\mp }\mathcal{A}$, respectively. Provided that the propagation constant is $k=1$, and ${\lambda }_{-}=1$, ${\lambda }_{+}=2$, ${\gamma }_{-}=0$, ${\gamma }_{+}=1$ are the weighting coefficients in the Robin condition, the relative boundary error $e_N$ exhibits the behavior shown in Figure 5. As it appears from Figure 6, the selection of the expansion order $N=10$ results in an accurate Fourier-Hankel series representation $v_N(x,y)$ of the solution (with error $e_N<1\mathrm{\%}$). The spatial distribution of $v_N(x,y)$ is shown in Figure 7, whereas the logarithmic magnitude of the relevant expansion coefficients $A_{p,m}$ and $B_{p,m}$ ($p=1,2$) is plotted in Figure 8.

Figure 5

Relative boundary error ${\mathit{e}}_{\mathit{N}}$ as a function of order N of the truncated Fourier-Hankel series expansion representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{1}\mathbf{/}\mathbf{2}$ , ${\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{-}}\mathbf{=}\mathbf{3}$ , ${\mathit{k}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}^{\mathbf{+}}\mathbf{=}\mathbf{5}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{9}\mathbf{/}\mathbf{4}$ , and ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{2}\mathbf{(}\mathbf{1}\mathbf{-}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}{\mathbf{log}}_{\mathbf{2}}\mathbf{cos}\frac{\mathit{\pi }}{{\mathit{k}}^{\mathbf{+}}}\mathbf{)}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}\mathbf{\to }\mathbf{+}\mathbf{\infty }$ .

Figure 6

Boundary behavior along ${\mathit{\partial }}_{\mathbf{-}}\mathcal{A}$ (a) and ${\mathit{\partial }}_{\mathbf{+}}\mathcal{A}$ (b) of the partial sum ${\mathit{U}}_{\mathit{N}}\mathbf{(}\mathit{\varrho }\mathbf{,}\mathit{\vartheta }\mathbf{)}$ of order $\mathit{N}\mathbf{=}\mathbf{10}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{1}\mathbf{/}\mathbf{2}$ , ${\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{-}}\mathbf{=}\mathbf{3}$ , ${\mathit{k}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}^{\mathbf{+}}\mathbf{=}\mathbf{5}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{9}\mathbf{/}\mathbf{4}$ , and ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{2}\mathbf{(}\mathbf{1}\mathbf{-}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}{\mathbf{log}}_{\mathbf{2}}\mathbf{cos}\frac{\mathit{\pi }}{{\mathit{k}}^{\mathbf{+}}}\mathbf{)}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}\mathbf{\to }\mathbf{+}\mathbf{\infty }$ .

Figure 7

Spatial distribution of the Fourier-Hankel series expansion ${\mathit{v}}_{\mathit{N}}\mathbf{(}\mathit{x}\mathbf{,}\mathit{y}\mathbf{)}$ of order $\mathit{N}\mathbf{=}\mathbf{10}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{1}\mathbf{/}\mathbf{2}$ , ${\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{-}}\mathbf{=}\mathbf{3}$ , ${\mathit{k}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}^{\mathbf{+}}\mathbf{=}\mathbf{5}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{9}\mathbf{/}\mathbf{4}$ , and ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{2}\mathbf{(}\mathbf{1}\mathbf{-}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}{\mathbf{log}}_{\mathbf{2}}\mathbf{cos}\frac{\mathit{\pi }}{{\mathit{k}}^{\mathbf{+}}}\mathbf{)}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}\mathbf{\to }\mathbf{+}\mathbf{\infty }$ .

Figure 8

Logarithmic magnitude of the coefficients ${\mathit{A}}_{\mathit{p}\mathbf{,}\mathit{m}}$ and ${\mathit{B}}_{\mathit{p}\mathbf{,}\mathit{m}}$ ( $\mathit{p}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{2}$ ) relevant to the Fourier-Hankel series expansion ${\mathit{v}}_{\mathit{N}}\mathbf{(}\mathit{x}\mathbf{,}\mathit{y}\mathbf{)}$ representing the solution of the Robin problem for the Helmholtz equation in the supershaped annulus described by the Gielis formula with parameters ${\mathit{k}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{1}\mathbf{/}\mathbf{2}$ , ${\mathit{d}}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}\mathbf{3}\mathbf{/}\mathbf{4}$ , ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{-}}\mathbf{=}\mathbf{2}$ , ${\mathit{\nu }}_{\mathit{y}}^{\mathbf{-}}\mathbf{=}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{-}}\mathbf{=}\mathbf{3}$ , ${\mathit{k}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}{\mathit{k}}^{\mathbf{+}}\mathbf{=}\mathbf{5}$ , ${\mathit{d}}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{d}}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{9}\mathbf{/}\mathbf{4}$ , and ${\mathit{\nu }}_{\mathit{x}}^{\mathbf{+}}\mathbf{=}{\mathit{\nu }}_{\mathit{y}}^{\mathbf{+}}\mathbf{=}\mathbf{2}\mathbf{(}\mathbf{1}\mathbf{-}{\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}{\mathbf{log}}_{\mathbf{2}}\mathbf{cos}\frac{\mathit{\pi }}{{\mathit{k}}^{\mathbf{+}}}\mathbf{)}$ , ${\mathit{\nu }}_{\mathbf{0}}^{\mathbf{+}}\mathbf{\to }\mathbf{+}\mathbf{\infty }$ .

Remark It has been observed that $L^2$ norm of the difference between the exact solution and the relevant approximation is generally negligible. Point-wise convergence seems to be verified in the considered domains, with the only exception of a set of measure zero consisting of quasi-cusped points. In the neighborhood of these points, oscillations of the truncated order solution, recalling the classical Gibbs phenomenon, usually take place (see Figure 6).

Conclusion

A Fourier-like projection method, in combination with the adoption of a suitable stretched coordinate system, has been developed for solving the Robin problem for the Helmholtz equation in supershaped annuli. In this way, analytically based expressions of the solution of the considered class of BVP s can be derived by using classical quadrature rules, so overcoming the need for cumbersome numerical techniques such as finite-difference or finite-element methods. The proposed approach has been successfully validated by means of a dedicated numerical procedure based on the computer-aided algebra tool Mathematica^©. A point-wise convergence of the expansion series representing the solution seems to be verified with the only exception of a set of measure zero consisting of the quasi-cusped points along the boundary of the problem domain. In these points, Gibbs-like oscillations may occur. The computed results are found to be in good agreement with the theoretical findings on a Fourier series.