Fokas method

Suppose that ${\displaystyle u}$ and ${\displaystyle v}$ both satisfy Laplace's equation in the interior of a convex bounded polygon ${\displaystyle \Omega }$. It follows that

${\displaystyle 0=v(u_{xx}+u_{yy})-u(v_{xx}+v_{yy})={\frac {\partial }{\partial x}}{\big (}vu_{x}-uv_{x}{\big )}+{\frac {\partial }{\partial y}}{\big (}vu_{y}-uv_{y}{\big )}.}$

Then Green's theorem implies the relation

${\displaystyle \int _{\partial \Omega }\left[\left(vu_{x}-uv_{x}\right)dy-\left(vu_{y}-uv_{y}\right)dx\right]=0.}$

(Eq.14)

In order to express the integrand of the above equation in terms of just the Dirichlet and Neumann boundary values, we parameterize ${\displaystyle u(x,y)}$ and ${\displaystyle v(x,y)}$ in terms of the arc length, ${\displaystyle s}$, of ${\displaystyle \partial \Omega }$. This leads to

${\displaystyle u_{x}dy-u_{y}dx={\frac {\partial u}{\partial w}}ds,}$

(Eq.15)

where ${\displaystyle {\frac {\partial u}{\partial w}}}$ denotes the normal derivative.

In order to further simplify the global relation, we introduce the complex variable ${\displaystyle z=x+iy}$, and its conjugate ${\displaystyle {\bar {z}}=x-iy}$. We then choose the test function ${\displaystyle v=e^{-i\lambda z}}$, leading to the global relation for Laplace's equation:

${\displaystyle \int _{\partial \Omega }e^{-i\lambda z}\left[{\frac {\partial u}{\partial w}}+\lambda u{\frac {dz}{ds}}\right]ds=0,\qquad \lambda \in \mathbb {C} .}$

(Eq.16)

A similar argument can also be used in the presence of a forcing term (giving a non-zero right-hand side). An identical argument works for the Helmholtz equation

${\displaystyle u_{xx}+u_{yy}+4\beta ^{2}u=0}$

and the modified Helmholtz equation

${\displaystyle u_{xx}+u_{yy}-4\beta ^{2}u=0.}$

Choosing respective test functions ${\displaystyle v=e^{-i\beta (\lambda z+{\frac {\bar {z}}{\lambda }})}}$ and ${\displaystyle v=e^{-i\beta (\lambda z-{\frac {\bar {z}}{\lambda }})}}$ lead to respective global relations

${\displaystyle \int _{\partial \Omega }e^{-i\beta (\lambda z+{\frac {\bar {z}}{\lambda }})}\left[{\frac {\partial u}{\partial w}}+\beta \left(\lambda {\frac {dz}{ds}}-{\frac {1}{\lambda }}{\frac {d{\bar {z}}}{ds}}\right)u\right]ds=0,\qquad \lambda \in \mathbb {C} \backslash \{0\},}$

and

${\displaystyle \int _{\partial \Omega }e^{-i\beta (\lambda z-{\frac {\bar {z}}{\lambda }})}\left[{\frac {\partial u}{\partial w}}+\beta \left(\lambda {\frac {dz}{ds}}+{\frac {1}{\lambda }}{\frac {d{\bar {z}}}{ds}}\right)u\right]ds=0,\qquad \lambda \in \mathbb {C} \backslash \{0\}.}$

These three cases deal with more general second order elliptic constant coefficient PDEs through a suitable linear change of variables.

The Dirichlet to Neumann map for a convex polygon Suppose that ${\displaystyle \Omega }$ is the interior of a bounded convex polygon specified by the corners ${\displaystyle z_{1},z_{2},\dots \ z_{n}}$. In this case, the global relation Eq.16 takes the form

${\displaystyle \sum _{j=1}^{n}u_{j}(\lambda )=0,\ \lambda \in \mathbb {C} ,}$

(Eq.17)

where

${\displaystyle u_{j}(\lambda )={\frac {1}{2}}\int _{z_{j}}^{z_{j}+1}e^{-i\lambda z}(u_{x}-iu_{y})\,dz,}$

(Eq.18)

or

${\displaystyle u_{j}(\lambda )=\int _{z_{j}}^{z_{j}+1}e^{-i\lambda z}\left({\frac {\partial u}{\partial w}}+\lambda u{\frac {dz}{ds}}\right)\,ds.}$

(Eq.19)

The side ${\displaystyle S_{j}}$, which is the side between ${\displaystyle z_{j}}$ and ${\displaystyle z_{j+1}}$, can be parametrized by

${\displaystyle z(t)=m_{j}+th_{j},\ m_{j}={\frac {z_{j}+z_{j+1}}{2}},\ h_{j}={\frac {z_{j+1}-z_{j}}{2}},\ t\in [-1,1].}$

Hence,

${\displaystyle ds=|h_{j}|\,dt,\ dz=h_{j}\,dt.}$

The functions ${\displaystyle u}$ and ${\displaystyle {\frac {\partial u}{\partial w}}}$ can be approximated in terms of Legendre polynomials:

${\displaystyle u^{j}(t)\approx \sum _{\ell =0}^{N-1}a_{\ell }^{j}P_{\ell }(t),\quad {\frac {\partial u^{j}(t)}{\partial w}}\approx \sum _{\ell =0}^{N-1}b_{\ell }^{j}P_{\ell }(t),}$

(Eq.20)

where for the cases of the Dirichlet, Neumann or Robin boundary value problems either ${\displaystyle a_{\ell }^{j}}$, ${\displaystyle b_{\ell }^{j}}$or a linear combination of ${\displaystyle a_{\ell }^{j}}$and ${\displaystyle b_{\ell }^{j}}$ is given.

Equation Eq.19 now becomes an approximate global relation, where

${\displaystyle u_{j}(\lambda )=\sum _{\ell =0}^{N-1}e^{-im_{j}\lambda }(b_{\ell }^{j}|h_{j}|+\lambda a_{\ell }^{j}h_{j}){\widehat {P}}_{\ell }(\lambda h_{j}),\ \lambda \in \mathbb {C} ,}$

(Eq.21)

with ${\displaystyle {\hat {P}}_{\ell }(\lambda )}$ denoting the Fourier transform of ${\displaystyle P_{\ell }(t)}$, i.e.,

${\displaystyle {\widehat {P}}_{\ell }(\lambda )=\int _{-1}^{1}e^{-i\lambda t}P_{\ell }(t)\,dt,\quad \lambda \in \mathbb {C} .}$

(Eq.22)

${\displaystyle {\hat {P}}_{\ell }(\lambda )}$ can be computed numerically via ${\displaystyle {\hat {P}}_{\ell }(\lambda )={\frac {\sqrt {-2i\pi \lambda }}{-i\lambda }}I_{\ell +{\frac {1}{2}}}(-i\lambda ),}$ where ${\displaystyle I_{\ell }}$ denotes the modified Bessel function of the first kind.

The global relation involves ${\displaystyle N}$ unknown constants (for the Dirichlet problem, these constants are ${\displaystyle b_{\ell }^{j}}$). By evaluating the global relation at a sufficiently large number of different values of ${\displaystyle \lambda }$, the unknown constants can be obtained via the solution of a system of algebraic equations.

It is convenient to choose the above values of ${\displaystyle \lambda }$ on the ${\displaystyle n}$ rays ${\displaystyle -{\bar {h}}_{j}\rho ,\ \rho >0,\ j=1,\dots ,n.}$ For this choice, as ${\displaystyle |\lambda |\rightarrow \infty }$, the relevant system is diagonally dominant, thus its condition number is very small.[3]

Whilst the global relation is valid for non-convex domains ${\displaystyle \Omega }$, the above collocation method becomes numerically unstable.[4] A heuristic explanation for this ill-conditioning in the case of the Laplace equation is as follows. The `test functions' ${\displaystyle e^{-i\lambda z}}$ grow/decay exponentially in certain directions of ${\displaystyle \lambda }$. When using a sufficiently large selection of complex ${\displaystyle \lambda }$-values, located in all directions from the origin, each side of a convex polygon will for many of these ${\displaystyle \lambda }$-values encounter larger test functions than do the remaining sides. This is exactly the same argument that motivates the `ray' choice of collocation points given by ${\displaystyle -{\overline {h}}_{j}\rho }$, which yield a diagonally dominant system. In contrast, for a non-convex polygon, boundary regions in indented regions will always be dominated by effects from other boundary parts, no matter the ${\displaystyle \lambda }$-value. This can easily be overcome by splitting up the domain into numerous convex regions (introducing fictitious boundaries) and matching the solution and normal derivative across these internal boundaries. Such splitting also allows the extension of the method to exterior/unbounded domains ${\displaystyle \Omega }$ (see below).

Let ${\displaystyle G=G(z,\zeta )}$ be the associated fundamental solution of the PDE satisfied by ${\displaystyle u}$. In the case of straight edges, Green's representation theorem leads to

${\displaystyle {\begin{aligned}u(x,y)&=\sum _{j=1}^{n}\int _{-1}^{1}{\Big [}G(z,\zeta (t)){\frac {\partial u}{\partial w}}|h_{j}|+iu{\Big (}{\frac {\partial G}{\partial \zeta (t)}}(z,\zeta (t))h_{j}-{\frac {\partial G}{\partial {\bar {\zeta }}(t)}}(z,\zeta (t)){\bar {h}}_{j}{\Big )}{\Big ]}dt\\&\approx \sum _{j=1}^{n}\sum _{l=0}^{N-1}\int _{-1}^{1}{\Big [}G(z,\zeta (t))b_{l}^{j}P_{l}(t)|h_{j}|+ia_{l}^{j}{\Big (}{\frac {\partial G}{\partial \zeta (t)}}(z,\zeta (t))h_{j}-{\frac {\partial G}{\partial {\bar {\zeta }}(t)}}(z,\zeta (t)){\bar {h}}_{j}{\Big )}P_{l}(t){\Big ]}dt.\end{aligned}}}$

(Eq.23)

Due to the orthogonality of the Legendre polynomials, for a given ${\displaystyle z=x+iy}$, the integrals in the above representation are Legendre expansion coefficients of certain analytic functions (written in terms of ${\displaystyle G}$). Hence the integrals can be computed rapidly (all at once) by expanding the functions in a Chebyshev basis (using the FFT) and then converting to a Legendre basis.[5] This can also be used to approximate the `smooth' part of the solution after adding global singular functions to take care of corner singularities.

The method can be extended to variable coefficient PDEs and curved boundaries in the following manner (see [6]). Suppose that ${\displaystyle \alpha (x,y)\in \mathbb {C} ^{2\times 2}}$ is a matrix valued function, ${\displaystyle \beta (x,y)\in \mathbb {C} ^{2}}$ a vector valued function and ${\displaystyle \gamma }$ a function (all sufficiently smooth) defined over ${\displaystyle \Omega }$. Consider the formal PDE in divergence form:

${\displaystyle \nabla \cdot (\alpha \nabla u)+\nabla \cdot (\beta u)+\gamma u=0.}$

(Eq.24)

Assume that the domain ${\displaystyle \Omega }$ is a bounded connected Lipschitz domain whose boundary consists of a finite number of vertices connected by ${\displaystyle C^{1}}$ arcs. Denote the corners of ${\displaystyle \Omega }$ in anticlockwise order as ${\displaystyle \{z_{j}\}_{1}^{n}}$ with the side ${\displaystyle \Gamma _{j}}$, joining ${\displaystyle z_{j}}$ to ${\displaystyle z_{j+1}}$. ${\displaystyle \Gamma _{j}}$ can be parametrised by

${\displaystyle [-1,1]\ni t\rightarrow (x_{j}(t),y_{j}(t))\in \mathbb {R} ^{2},}$

where we assume that the parametrisation is ${\displaystyle C^{1}}$.

The adjoint of equation Eq.24 is given by

${\displaystyle \nabla \cdot (\alpha ^{T}\nabla v)-\beta \cdot \nabla v+\gamma v=0.}$

(Eq.25)

The expression ${\displaystyle v\times }$Eq.24${\displaystyle -u\times }$Eq.25 can be written in the form

${\displaystyle \nabla \cdot [v(\alpha \nabla u)-u(\alpha ^{T}\nabla v)+\beta uv]=0.}$

(Eq.26)

Integrating across the domain and applying the divergence theorem we recover the global relation (${\displaystyle n}$ denotes the outward normal):

${\displaystyle \int _{\partial \Omega }u{\big [}(n\cdot \beta )v-n\cdot (\alpha ^{T}\nabla v){\big ]}+v{\big [}n\cdot (\alpha \nabla u){\big ]}ds=0.}$

(Eq.27)

Define ${\displaystyle l_{j}(t)={\sqrt {{\dot {x}}_{j}(t)^{2}+{\dot {y}}_{j}(t)^{2}}}}$ along the curve ${\displaystyle \Gamma _{j}}$ and assume that ${\displaystyle l_{j}(t)>0}$. Suppose that we have a one-parameter family of solutions of the adjoint equation, ${\displaystyle v(x,y;\lambda )=v_{\lambda }}$, for some ${\displaystyle \lambda \in {\mathcal {C}}}$, where ${\displaystyle {\mathcal {C}}}$ denotes the collocation set. Denoting the solution ${\displaystyle u}$ alongside ${\displaystyle \Gamma _{j}}$ by ${\displaystyle u^{j}}$, the unit outward normal by ${\displaystyle n_{j}}$ and analogously the oblique derivative by ${\displaystyle n_{j}\cdot (\alpha \nabla u^{j})}$, we define the following important transform:

${\displaystyle {\hat {\mu }}_{j}(\lambda ):=\int _{-1}^{1}{\big \{}u^{j}{\big [}(n_{j}\cdot \beta )v_{\lambda }-n_{j}\cdot (\alpha ^{T}\nabla v_{\lambda }){\big ]}+v_{\lambda }{\big [}n_{j}\cdot (\alpha \nabla u^{j}){\big ]}{\big \}}l_{j}(t)dt.}$

(Eq.28)

Using Eq.28 , the global relation Eq.27 becomes

${\displaystyle \sum _{j=1}^{n}{\hat {\mu }}_{j}(\lambda )=0.}$

(Eq.29)

For separable PDEs, a suitable one-parameter family of solutions ${\displaystyle v_{\lambda }}$ can be constructed. If we expand each ${\displaystyle u^{j}}$ and its derivative ${\displaystyle n_{j}\cdot (\alpha \nabla u^{j})}$ along the boundary ${\displaystyle \Gamma _{j}}$ in Legendre polynomials, then we cover a similar approximate global relation as before. To compute the integrals that form the approximate global relation, we can use the same trick as before - expanding the function integrated against Legendre polynomials in a Chebyshev series and then converting to a Legendre series. A major advantage of the method in this scenario is that it is a boundary-based method which does not need any knowledge of the corresponding Green's function. Hence, it is more applicable than boundary integral methods in the setting of variable coefficients.

A major advantage of the above collocation method is that the basis choice (Legendre polynomials in the above discussion) can be flexibly chosen to capture local properties of the solution along each boundary. This is useful when the solution has different scalings in different regions of ${\displaystyle \Omega }$, but is particularly useful for capturing singular behavior, for example, near sharp corners of ${\displaystyle \Omega }$.

We consider the acoustic scattering problem solved in [7] by the method. The solution ${\displaystyle u}$ satisfies Helmholtz equation in ${\displaystyle \mathbb {R} ^{2}\backslash \{x\in [-1,1],y=0\}}$ with frequency ${\displaystyle k_{0}=2\beta }$, along with the Sommerfeld radiation condition at infinity:

${\displaystyle \lim _{r\rightarrow \infty }r^{1/2}\left({\frac {\partial }{\partial r}}-ik_{0}\right)u(x,y)=0,}$

(Eq.30)

where ${\displaystyle r={\sqrt {x^{2}+y^{2}}}}$. The boundary condition along the plate is

${\displaystyle {\frac {\partial u}{\partial y}}(x,0)=-{\frac {\partial u_{I}}{\partial y}}(x,0),\quad x\in [-1,1]}$

(Eq.31)

for the incident field

${\displaystyle u_{I}(x,y)=e^{-ik_{0}\cos \theta x-ik_{0}\sin \theta y}.}$

(Eq.32)

By considering the domains ${\displaystyle y>0}$ and ${\displaystyle y<0}$ separately and matching the global relations, the global relation for this problem becomes

${\displaystyle {\begin{aligned}&\int _{\mathbb {R} \backslash [-1,1]}e^{-i\beta x(\lambda +{\frac {1}{\lambda }})}u_{y}(x,0)dx\\&+\int _{[-1,1]}e^{-i\beta x(\lambda +{\frac {1}{\lambda }})}\left[u_{y}(x,0)+{\frac {\beta }{2}}\left(\lambda -{\frac {1}{\lambda }}\right)[u](x,0)\right]dx=0,\qquad \lambda \in \Lambda ,\end{aligned}}}$

(Eq.33)

with ${\displaystyle \Lambda =(-1,0)\cup (1,\infty )\cup \{e^{-i\theta }:0<\theta <\pi \}}$ and where ${\displaystyle [u](x,0)=u(x,0_{+})-u(x,0_{-})}$ denotes the jump in ${\displaystyle u}$ across the plate. The complex collocation points are allowed precisely because of the radiation condition. To capture the endpoint singularities, we expand ${\displaystyle [u](x,0)}$ for ${\displaystyle x\in [-1,1]}$ in terms of weighted Chebyshev polynomials of the second kind:

${\displaystyle C_{m}(x)={\sqrt {1-x^{2}}}U_{m}(x).}$

(Eq.34)

These have the following Fourier transform:

${\displaystyle \int _{-1}^{1}e^{i\lambda t}{\sqrt {1-t^{2}}}U_{m}(t)dt={\frac {(m+1)i^{m}\pi }{\lambda }}J_{m+1}(\lambda ),}$

(Eq.35)

where ${\displaystyle J_{\alpha }(\cdot )}$ denotes the Bessel function of the first kind of order ${\displaystyle \alpha }$. For the derivative ${\displaystyle u_{y}}$ along ${\displaystyle \mathbb {R} \backslash [-1,1]}$, a suitable basis choice are Bessel functions of fractional order (to capture the singularity and algebraic decay at infinity).

We introduce the dimensionless frequency ${\displaystyle {\tilde {k}}_{0}=k_{0}L^{-1}}$, where ${\displaystyle L}$ is the length of the plate. The figure below shows the convergence of the method for various ${\displaystyle {\tilde {k}}_{0}}$. Here ${\displaystyle N}$ is the number of basis functions ${\displaystyle C_{m}}$ used to approximate the jump ${\displaystyle [u]}$ across the plate. The maximum relative absolute error is the maximum error of the computed solution divided by the maximum absolute value of the solution. The figure is for ${\displaystyle \theta =\pi /3}$ and shows the quadratic-exponential convergence of the method, namely the error decreases like ${\displaystyle {\mathcal {O}}(\rho ^{N^{2}})}$ for some positive ${\displaystyle \rho <1}$. More complicated geometries (including different angles of touching boundaries and infinite wedges) can also be dealt with in a similar fashion as well as more complicated boundary conditions such as those modeling elasticity.[8][9]

Convergence results for the method and different ${\displaystyle {\tilde {k}}_{0}}$.