Spectral Method on Disk

REFERENCES

1. John P. Boyd, Fu Yu, Comparing seven spectral methods for interpolation and for solving the Poisson equation in a disk: Zernike polynomials, Logan–Shepp ridge polynomials, Chebyshev–Fourier Series, cylindrical Robert functions, Bessel–Fourier expansions, square-to-disk conformal mapping and radial basis functions. JCP(2011). https://doi.org/10.1016/j.jcp.2010.11.011
2. Wilber, H., Townsend, A., & Wright, G. B. Computing with functions in spherical and polar geometries II. The disk. SISC(2017). https://epubs.siam.org/doi/10.1137/16M1070207

谱方法在求解高维问题时面临着维数灾难的问题，即随着维数的增加，计算量和存储需求呈指数级增长，这迫使我们寻找更高效的高维区域的谱方法。最简单的高维区域是平面上的圆盘，这里考察一些常用的谱方法在圆盘区域上的表现。

为了逼近定义在单位圆盘上的函数，通常有以下几大类方法[2]：

1. 径向基函数(Radial Basis Functions, RBFs)：这类方法使用径向基函数来构建逼近，常见的RBFs包括高斯函数、多项式函数等，优势是能够处理不规则分布的数据点，适用于高维空间中的插值和逼近问题，但是对于这里考虑的圆盘问题，RBFs的精度受限于其较大的条件数，当求解中等规模的问题时将会丢失3到5位精度。
2. 共形变换(Conformal Mapping)：通过将圆盘区域映射到一个更简单的区域（如方形），在该区域上应用标准的谱方法，然后将结果映射回圆盘。这种方法的优点是可以利用现有的谱方法，但缺点是共形变换可能会引入数值不稳定性，尤其是在区域的角点处会出现新的奇异性(corner singularity)。
3. 基函数展开(Basis Expansions)：在径向和角向上使用不同的基函数进行展开，通常在角向上使用傅里叶级数，在径向上使用适合圆盘边界条件的函数（如Zernike多项式、Bessel函数等）。这种方法的优点是能够很好地适应圆盘的几何特性，但需要选择合适的基函数来满足边界条件和解析性要求。

这篇文章聚焦于基函数展开方法，比较了几种常用的基函数展开方法在圆盘区域上的表现，并分析了它们的优缺点。

对于单位圆盘上的问题，通常使用极坐标来参数化，下面我们使用以下形式来表示上的函数：

其中径向参数。这样的参数化对自身以及系数函数的性质提出了一些要求。首先，由于函数在方向上是周期性的，上述表示的圆盘函数必须满足

另一方面，因为笛卡尔坐标与极坐标之间满足与，所以频率为的项具有以下形式

由于函数在圆盘在原点处解析，而显然是全纯的，因此的解析性要求必须也在处解析，故存在非奇异函数和使得

又因为有，因此与必须满足与，所以

即与只包含偶数幂的项。由于性质与性质是解析函数内蕴的自然要求，因此在选择基函数进行展开时，我们希望所选的基函数能够满足这两个性质所要求的，以尽量保证逼近的函数具有正确的解析性。

Zernike Polynomials

Zernike多项式是定义在单位圆盘上的一类正交多项式，又被称为单侧Jacobi多项式，常用于光学和图像处理领域。它们可以表示为以下形式：

其中且，也即。被称为径向Zernike多项式，具体定义为

从定义式可以直接看出自然满足条件，因此适合用来逼近和。由于使用Jacobi多项式进行定义，该径向多项式自动继承了Jacobi多项式的性质，满足如下递推关系：

初始多项式为与。另外，在权重下正交，即

进一步有

其中时，当时。使用这一性质可知圆盘上的函数的Zernike展开

中的展开系数由以下二重积分给出：

实际计算时，可以使用数值积分方法来近似计算这些积分，在角向上可以使用FFT快速计算，在径向上通常使用Gauss积分来计算，即将上述积分转化为

再分别使用FFT和Gauss积分来数值计算内外两个积分。具体地来说，首先定义

其中或，对于一个固定的可以在复杂度下计算出的所有，而后续的径向积分变为

此处与分别为上述分段定义式中的常数系数与径向函数。使用Gauss积分需要个Gauss积分点与对应的权重，从而将上述积分近似为

因此计算出且所需的总复杂度为。当与具有相同数量级时，整体复杂度为，这对于中等规模的问题来说是可行的，但是不太适用于较大规模的问题。目前Zernike多项式的应用主要被缺少高效的快速算法所限制，虽然有一些文献提出了基于蝶形算法与HSS的快速算法，但通常具有较大的预处理复杂度，一般适用于时间发展问题这类需要频繁重复计算上述积分的情形，在本文考虑的问题中仍然不具备较大优势。上述计算的瓶颈在于径向上缺少类似于FFT的快速算法。

在实际计算中，通常需要将展开式进行截断，最常使用的是三角截断，即仅保留且的项，得到如下的近似表达式：

上述表达式中使用到的Zernike多项式的数量为，可以排列为以下的三角形

下图展示了不同阶数的Zernike多项式在圆盘上的分布情况。随径向阶数 增加(从上往下)可以看到径向波动增加，即沿着垂直方向从上往下看，图案变得越来越复杂，内部的“同心圆环”数量在增加。随角向频率增加(从中心向两侧)图像的角向“花瓣”增加。在同一行中，关于中心轴对称的左右两个图案(相同，互为相反数)的形状和色彩分布规律完全相同，仅仅是发生了一个固定角度的旋转。沿着金字塔最左侧和最右侧的边缘(满足的多项式)处的图案只包含从中心辐射出的“花瓣”，没有任何内部的同心圆环，此时的径向多项式退化为简单的幂函数，不再有内部的径向零点，波前畸变完全由角向的三角函数主导。

前21个Zernike多项式在圆盘上的分布情况

Poisson方程的Zernike谱方法

下面使用Zernike多项式作为基函数求解圆盘区域上的Poisson方程

其中是定义在上的一个函数，要求它具有下述形式的Zernike展开

现在我们希望找到一个Zernike展开形式的近似解

首先将Laplace算子表示为极坐标下的形式

注意到这一算子的可分离性，我们可以将上述方程中的与的展开式代入，得到

为了使上式成立，我们要求每一个频率对应的内层求和的项都必须相等，即

此时原本的二维问题被解耦为了一系列一维ODE

现在可以利用径向Zernike多项式的性质来为这些ODE构造稀疏求解器。