Non-singular Green's functions for the unbounded Poisson equation in one, two and three dimensions

TL;DR

This paper derives non-singular Green's functions for the unbounded Poisson equation in 1D, 2D, and 3D using a Fourier domain cutoff, which are useful for applications with a minimum length scale.

Contribution

It introduces a novel method to obtain non-singular Green's functions for the Poisson equation across multiple dimensions using a Fourier cutoff approach.

Findings

Derived explicit non-singular Green's functions for 1D, 2D, and 3D.

Provided the gradient of the non-singular Green's functions for vector field applications.

Applicable to problems with a minimum resolved length scale, such as mesh-based simulations.

Abstract

This paper is a revised version of the original paper of same title--published in Applied Mathematics Letters 89--containing some corrections and clarifications to the original text. We derive non-singular Green's functions for the unbounded Poisson equation in one, two and three dimensions, using a cut-off function in the Fourier domain to impose a smallest length scale when deriving the Green's function. The resulting non-singular Green's functions are relevant to applications which are restricted to a minimum resolved length scale (e.g. a mesh size h) and thus cannot handle the singular Green's function of the continuous Poisson equation. We furthermore derive the gradient vector of the non-singular Green's function, as this is useful in applications where the Poisson equation represents potential functions of a vector field.

Equations

$${\nabla }^{2}A(\mathbf{x})=-B(\mathbf{x})\text{where}A(\mathbf{x})\to 0\text{for}|\mathbf{x}|\to \infty$$

$$\mathbf{v}(\mathbf{x})=\nabla \times \mathbf{\psi }(\mathbf{x})-\nabla \varphi (\mathbf{x})\text{mbox}where\nabla \cdot \mathbf{\psi }(\mathbf{x})=0$$

$$\vartheta (\mathbf{x})=-{\nabla }^2\varphi (\mathbf{x})\text{mbox}and\mathbf{\omega }(\mathbf{x})=-{\nabla }^2\mathbf{\psi }(\mathbf{x})$$

$$\begin{array}[]{l c l}-k^{2}\widehat{\phi}(\boldsymbol{k})=-\widehat{\vartheta}(\boldsymbol{k})&\Leftrightarrow&\widehat{\phi}(\boldsymbol{k})=\widehat{G}(\boldsymbol{k})\,\widehat{\vartheta}(\boldsymbol{k})\\[8.5359pt] -k^{2}\widehat{\boldsymbol{\psi}}(\boldsymbol{k})=-\widehat{\boldsymbol{\omega}}(\boldsymbol{k})&\Leftrightarrow&\widehat{\boldsymbol{\psi}}(\boldsymbol{k})=\widehat{G}(\boldsymbol{k})\,\widehat{\boldsymbol{\omega}}(\boldsymbol{k})\end{array}\quad\text{where}\quad\widehat{G}(\boldsymbol{k})=\frac{1}{k^{2}}\quad\text{for}\quad k>0$$

$$\begin{array}[]{l c l}-k^{2}\,\widehat{\boldsymbol{v}}_{\vartheta}(\boldsymbol{k})=-\text{i}\,\boldsymbol{k}\,\widehat{\vartheta}(\boldsymbol{k})&\Leftrightarrow&\widehat{\boldsymbol{v}}_{\vartheta}(\boldsymbol{k})=\widehat{\boldsymbol{K}}(\boldsymbol{k})\,\widehat{\vartheta}(\boldsymbol{k})\\[8.5359pt] -k^{2}\,\widehat{\boldsymbol{v}}_{\omega}(\boldsymbol{k})=-\text{i}\,\boldsymbol{k}\times\widehat{\boldsymbol{\omega}}(\boldsymbol{k})&\Leftrightarrow&\widehat{\boldsymbol{v}}_{\omega}(\boldsymbol{k})=\widehat{\boldsymbol{K}}(\boldsymbol{k})\times\widehat{\boldsymbol{\omega}}(\boldsymbol{k})\end{array}\quad\text{where}\quad\widehat{\boldsymbol{K}}(\boldsymbol{k})=\frac{\text{i}\,\boldsymbol{k}}{k^{2}}\quad\text{for}\quad k>0$$

$$\widehat{\zeta}(k)=\left\{\begin{array}[]{l l l}1&\text{for}&0<k\leq\dfrac{\pi}{h}\\[5.69046pt] 0&\text{for}&k>\dfrac{\pi}{h}\end{array}\right.$$

$$\hat{G}(k)=\frac{\hat{\zeta }(k)}{k^2}\iff {\nabla }^{2}G(r)=-\zeta (r)$$

$$\hat{\zeta }(s)=(2\pi {\sigma }^2{)}^{\frac{d}{2}}{\int }_0^{\infty }\zeta (\rho )\frac{J_{\frac{d}{2}-1}(s\rho )}{(s\rho {)}^{\frac{d}{2}-1}}{\rho }^{d-1}d\rho$$

$$\widehat{\zeta}(s)=\displaystyle\int_{0}^{\infty}\frac{J_{\frac{d}{2}}(\rho)J_{\frac{d}{2}-1}(s\rho)}{s^{\frac{d}{2}-1}}\,d\rho=\left\{\begin{array}[]{l l l}1&\text{for}&0<s<1\\[5.69046pt] \frac{1}{2}&\text{for}&s=1\\[5.69046pt] 0&\text{for}&1<s\end{array}\right.$$

$$\sigma \frac{\pi }{h}=1\iff \sigma =\frac{h}{\pi }\Rightarrow s=\frac{hk}{\pi }\text{and}\rho =\frac{\pi r}{h}$$

$$\zeta (\rho )=\frac{J_{\frac{d}{2}}(\rho )}{(2\pi {\sigma }^2\rho {)}^{\frac{d}{2}}}$$

$$G(\rho)=-\dfrac{\sigma}{\pi}\big{(}\text{Si}(\rho)\,\rho+\cos(\rho)\big{)}+C_{1}$$

$$C_1=\frac{1}{2}L\iff G(\rho )\to -\frac{1}{2}(r-L)\text{for}\rho \to \infty$$

$$K(\rho )=-\frac{\text{Si}(\rho )}{\pi }$$

$$G(\rho )=-\frac{1}{2\pi }\text{Bi}(\rho )+C_2$$

$$C_2=-\frac{1}{2\pi }\left(\text{ln}\left(\frac{2\sigma }{L}\right)-\gamma \right)\iff G(\rho )\to -\frac{1}{2\pi }\text{ln}\left(\frac{r}{L}\right)\text{for}\rho \to \infty$$

$$\mathbf{K}(\rho )=-K(\rho ){\mathbf{e}}_r\text{with}K(\rho )=-\frac{1}{\sigma }\frac{dG}{d\rho }=\frac{1}{2\pi \sigma }\frac{1-J_0(\rho )}{\rho }$$

$$G(\rho )=\frac{1}{2{\pi }^2\sigma }\frac{\text{Si}(\rho )}{\rho }+C_3\text{with}G(0)=\frac{1}{2{\pi }^2\sigma }$$

$$C_3=0\iff G(\rho )\to \frac{1}{4\pi r}\text{for}\rho \to \infty$$

$$\mathbf{K}(\rho )=-K(\rho ){\mathbf{e}}_r\text{with}K(\rho )=-\frac{1}{\sigma }\frac{dG}{d\rho }=\frac{1}{2{\pi }^2{\sigma }^2}\frac{\text{Si}(\rho )-\sin (\rho )}{{\rho }^2}$$

Full text

We derive non-singular Green’s functions for the unbounded Poisson equation in one, two and three dimensions, using a cut-off function in the Fourier domain to impose a smallest length scale when deriving the Green’s function. The resulting non-singular Green’s functions are relevant to applications which are restricted to a minimum resolved length scale (e.g. a mesh size $h$) and thus cannot handle the singular Green’s function of the continuous Poisson equation. We furthermore derive the gradient vector of the non-singular Green’s function, as this is useful in applications where the Poisson equation represents potential functions of a vector field.

††journal: arXiv

1 Introduction

The use of Green’s functions for solving linear differential equations in an unbounded domain, i.e. with free-space boundary conditions, is a frequently used methodology. The Green’s function represents the impulse response function to a linear differential operator, such as the Laplace operator in the case of the Poisson equation. Once obtained, the Green’s function can be used, by utilizing the superposition principle to obtain the solution of an inhomogeneous equation by convolving the Green’s function with the right-hand-side field of the equation.

The Green’s function of the Laplace operator is singular at its origin. Applied in discretized numerical calculations, the singularity of the Green’s function evidently causes a number of difficulties. In order to amend this, smoothing regularization techniques have been applied (e.g. [2, 3, 4]) which imposes the effect of a continuous and smooth field distribution around the discrete points, which in effect can be used to avoid the singularity of the Green’s function. However, most regularization methods that have been applied are based on functions that only conserve a finite number of field moments, and are thus only accurate up to a finite order of convergence rate.

Vico et al. [5] derived non-singular Green’s function in the two and three dimensional Fourier domain by imposing an isotropic maximum length scale of the real domain. This method was shown to provide a spectral accuracy in computations, where the real domain was sufficiently extended in all directions to provide a sufficient resolution when evaluating the Green’s function in the Fourier domain. However, as the Green’s function is a radial function the imposed maximum length scale dictates the resolution criteria in the Fourier domain, and thus the length of the integrated domain, for all directions. This makes the method of Vico et al. [5] potentially inefficient for elongated domains, as the imposed length scale must be proportional to the maximum length of the domain to insure that the integration includes the full domain.

In this work we impose a minimum length scale in the real domain, analogous to a discretization length, in order to obtain a non-singular Green’s functions in the real domain. This is done using a cut-off function in the Fourier domain to impose a minimum length scale, allowing us to derive the Green’s function analytically by Fourier analysis. As no explicit smoothing function is applied, the obtained non-singular Green’s function gives an optimal accuracy when used in numerical calculations, as it is only subject to quadrature errors, converging in a rate proportional to the smoothness of the right-hand-side field of the Poisson equation.

As with the method of Vico et al. [5] the presented method also imposes an equivalent length scale for all directions, albeit here as a minimum length scale of the solution. Consequently, if the discretization is non-isotropic the proposed method will provide a solution, in all directions, which limited to the minimum length scale of the lowest resolved direction.

The gradient vector of the non-singular Green’s function is also derived. This function is useful in applications where the Poisson equation represents a vector or scalar potential function of a vector field as it enables a direct solution of the vector field.

The proposed non-singular Green’s functions may be used directly in efficient computational methods such as the fast multipole method [6, 8] or the mesh based FFT solver [7]. The mesh based Poisson solver was recently shown capable of obtaining a high order convergence rate by using regularized Green’s functions [9, 10, 11]. The mesh based solver was extended to handle domains with a combination of unbounded and periodic directions by Chatelain and Koumoutsakos [12] and for high order accuracy in Spietz et al. [13]. Here it is utilized that the domain is intrinsically bounded in the periodic directions, yielding a straightforward Green’s function in the Fourier domain, which may be combined with the Fourier transformed unbounded Green’s function of reduced dimensionality for the unbounded directions.

Using the non-singular Green’s functions presented here in the mesh based Poisson solver effectively results in a solution of spectral accuracy for both unbounded domains as well as a domains with a combination of unbounded and periodic directions. The reader is referred to the open source software [14] for a numerical implementation of this.

2 Methodology

The Poisson equation in an unbounded domain is formally stated as:

[TABLE]

Here $B(\boldsymbol{x})$ is a known compactly supported field and $A(\boldsymbol{x})$ is the desired solution field. In many applications such as astrophysics, electrodynamics and vortex dynamics, the vector field to be solved $\boldsymbol{v}(\boldsymbol{x})$ is described by potential functions using the Helmholtz decomposition:

[TABLE]

The fundamental operations describing the pointwise conservation of the flux and circulation of the vector field $\boldsymbol{v}$ is the divergence $\vartheta(\boldsymbol{x})=\boldsymbol{\nabla}\cdot\boldsymbol{v}(\boldsymbol{x})$ and the curl $\boldsymbol{\omega}(\boldsymbol{x})=\boldsymbol{\nabla}\times\boldsymbol{v}(\boldsymbol{x})$, respectively. From Eq. (2) it follows that these may be expressed by the potential functions as:

[TABLE]

Both potentials may thus be obtained by solving a Poisson equation. By utilizing the linearity of the Poisson equation, and considering the equations transformed into the Fourier domain, we obtain the algebraic equations:

[TABLE]

Here  $\widehat{\cdot}$  denotes the field generated by the Fourier transform, $\boldsymbol{k}$ is the angular wave-number vector of the Fourier domain, and $k=|\boldsymbol{k}|$ is its radial length. $\widehat{G}(\boldsymbol{k})$ is the Green’s function in the Fourier domain representation. We may furthermore obtain the vector field of Eq. (2) directly by incorporating the gradient operator into the Green’s function as $\boldsymbol{K}(\boldsymbol{x})=\boldsymbol{\nabla}G(\boldsymbol{x})$, after which we obtain:

[TABLE]

As the Green’s function in the Fourier domain (Eq. (4)) does not have a compact support, an unbounded Fourier domain is needed to obtain the exact solution. However, for a discrete approximation, such as that used in numerical simulations, the discretization of a function in the real domain, imposes an upper bound $k_{s}$ of the function image in the Fourier domain. This upper bound is determined by the Nyquist-Shannon sampling theorem as $k_{s}=\pi/h$ being the highest resolved angular wave-number corresponding to the discretization length $h$. Thus, the exact Green’s function of Eq. (4) cannot not be represented in a discretized domain, and using it will introduce errors in the calculation resulting in a limited convergence rate of $\mathcal{O}(h^{2})$ [15]. In order to amend this and derive a Green’s function which is bounded by a maximum wavenumber, we regularize the Green’s function using a radial cut-off function in the Fourier domain:

[TABLE]

This leads to a regularized Green’s function and a corresponding equation for the impulse response function:

[TABLE]

where $r=|\boldsymbol{x}|$ is the radial length of the position vector. Using the radial cut-off function of Eq. (6), it is seen that the Green’s function of the non-regularized equation (Eq. (4)) is unchanged for $0<k\leq\pi/h$.

3 Results

By utilizing the radial symmetry of the Green’s function, the regularization function in the real domain $\zeta(\rho)$ may be obtained by considering the $d$-dimensional radial Fourier transform [16] for $d\geq 1$:

[TABLE]

where we have introduced the normalized coordinates $s=\sigma k$ and $\rho=r/\sigma$ relative to the length scale $\sigma$. The function $J_{\nu}$ is the Bessel function of the first kind and of order $\nu$. Using an integral property of the Bessel functions of the first kind [17], we may satisfy Eq. (6) by:

[TABLE]

The scaling length $\sigma$ is now determined by fulfilling $s=1$ for $k=k_{s}=\pi/h$:

[TABLE]

Here $h$ is the discretization length, chosen so the smallest physical length scale of the problem is sufficiently resolved.

Considering Eqs. (8) and (9), we directly see that the regularization function is given by:

[TABLE]

For the case of $d=2$ this is identical to a regularization function presented in the references [2, 4].

We may now obtain the regularized Green’s function in the real domain by radial integration of Eq. (7) in the unbounded real domain. For the one-dimensional case ($d=1$), we obtain:

[TABLE]

Here $\text{Si}(\rho)=\int_{0}^{\rho}\frac{\sin(t)}{t}\ dt$ is the sine integral function and $C_{1}$ is an integration constant, which we may use to define a reference value. In this work we determine $C_{1}$ such that we obtain an asymptotic behavior towards the non-regularized Green’s function for large $\rho$ whatever the value of $\sigma$:

[TABLE]

Here $L$ is an arbitrary reference length which typically represents the largest physical length scale of the problem, by which $L\gg h$.

The gradient vector of the Green’s function, which may be used to combine the gradient operator directly in Eq. (5), is in the one-directional case ($d=1$) given by:

[TABLE]

For the two-dimensional case ($d=2$) the radial integration yields the real domain Green’s function:

[TABLE]

where $\text{Bi}(\rho)=\int_{0}^{\rho}\frac{1-J_{0}(t)}{t}\ dt$ and $C_{2}$ is an integration constant, which we may use to define a reference value. The integral function Bi can be expressed by $\text{Bi}(\rho)=\text{Ji}_{0}(\rho)+\ln\left(\frac{\rho}{2}\right)+\gamma$ with $\gamma$ being the Euler constant and $\text{Ji}_{0}(\rho)=\int_{\rho}^{\infty}\frac{J_{0}(t)}{t}\ dt$ denoting the Bessel integral function which for a numerical calculation may be approximated efficiently by Chebyshev polynomials (see e.g. [18]). Determining the integration constant $C_{2}$ such that we obtain an asymptotic behavior towards the singular Green’s function for large $\rho$ whatever the value of $\sigma$:

[TABLE]

The gradient vector of the Green’s function, which may be used to calculate the vector field directly in Eq. (5), is in the two-directional case ($d=2$) given by:

[TABLE]

Here $K(\rho)\geq 0$ with $K(0)=0$ and $\boldsymbol{e}_{r}=\boldsymbol{x}/r$ is the radial unit vector. This result is equivalent to that presented in the references [2, 4].

Using the same approach for the three-dimensional domain ($d=3$), we obtain the Green’s function:

[TABLE]

The integration constant $C_{3}$ is again chosen such that we obtain an asymptotic behavior towards the singular Green’s function for large $\rho$ whatever the value of $\sigma$:

[TABLE]

The gradient vector of the Green’s function is in the three-directional case ($d=3$) obtained as:

[TABLE]

where we also have that $K(\rho)\geq 0$ and $K(0)=0$.

The derived non-singular Green’s functions are compared to the singular Green’s functions in Fig. 1. The oscillatory behavior of the non-singular Green’s functions is analogous to the sinc function approximation of the Dirac delta function, and ensures the conservation of the field moments.

For validation we have implementing the non-singular Green’s functions into the numerical methodology used in the Poisson solver presented in references [9, 10, 11, 13]. Using this method we are able to obtain an accuracy in the solution limited only by the machine precision, when tested on the same benchmark cases as is used in [9] (not shown). It is here emphasized that this order of precision is only feasible given a sufficiently smooth and resolved right-hand-side field of the Poisson equation. For less smooth and resolved fields, the presented non-singular Green’s functions still provides the highest possible accuracy allowed by the quadrature error. For the actual numerical implementation that we have used for these tests, the reader is referred to the provided open source software [14].

4 Conclusion

We have derived the non-singular Green’s functions for the unbounded Poisson equation, in one, two and in three dimensions, using a Fourier analysis. By applying a radial cut-off function in the Fourier domain, the regularized Green’s function is derived analytically in unbounded real domains subject to a minimum resolved length scale (e.g. a mesh size $h$). As a sharp cut-off function is used, smoothing errors are avoided and the obtained regularized Green’s function achieves an optimal accuracy when used to solve the Poisson equation on discretized fields.

The gradient vector of the Green’s function was furthermore derived. The obtained Green’s function for the two-dimensional case corresponds to the one stated in references [2, 4]. In the present work, we formally showed the derivation of this function by Fourier analysis, and we extended it to the one and three-dimensional cases.

Acknowledgement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. We would like to acknowledge the helpful support of Henrik Juul Spietz in developing the numerical methodology used in the referenced mesh based Poisson solver, and Denis-Gabriel Caprace for pointing out some needed improvements to the original manuscript.