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构造稀疏求解器，最终将得到一个带状线性系统。

为了简化记号，下面固定一个角向频率，并记。该频率对应的径向微分算子为

假设右端项的Zernike展开截断到总次数。令以及，则频率对应的径向右端项为

由于所有径向Zernike多项式都满足，若直接使用作为试探基，则齐次Dirichlet边界条件会转化为所有展开系数之和为零。为了将边界条件直接嵌入试探空间，定义

显然，而因子保证了该基函数在原点处具有频率所要求的正则性。

利用Jacobi多项式的参数转换公式，可以得到

其中。因此，每个都是两个相邻径向Zernike多项式之差，并且其最高次数为。

现在将频率对应的近似解写为

虽然右端项的最高次数不超过，但还包含次数为的Zernike多项式，因此近似解的最高次数可能达到。完整的二维近似解为。

对径向方程使用Galerkin方法，并取作为测试函数，得到

由于试探函数和测试函数在处均为零，并且在处满足相应的正则性条件，对左端分部积分后边界项消失，从而有

上述Dirichlet型Zernike基关于能量内积正交，并满足

因此，Laplace算子在该Galerkin基下对应一个对角刚度矩阵。

另一方面，由径向Zernike多项式的正交关系

可得

其中在截断末端约定。

定义未知向量和右端向量

令，并定义上双对角矩阵，其非零元素为

这里第一个公式适用于，第二个公式适用于。显式地，

于是，对于每一个角向频率，最终的Galerkin线性系统为

其中是对角矩阵，而是上双对角矩阵，因此每个径向系统都可以在线性复杂度下直接求解。由于，Galerkin系数还可以直接写为

若希望将结果重新表示为标准Zernike展开，则相应的系数满足

定义，则上述基变换可以写为，其中

矩阵的大小为，并且具有下双对角结构。因此标准Zernike系数满足

记，则是一个的三对角带状矩阵，它表示满足原点正则性条件和齐次Dirichlet边界条件的径向逆Laplace算子。

这一三对角结构也可以由逆Laplace算子的显式作用直接看出。对于最低阶径向模态，有

对于且，有

式和式右端各项的系数之和均为零，因此所得函数在处自动满足齐次Dirichlet边界条件。

具体而言，对于，右端系数只会对解系数、和产生贡献，相应的矩阵元素为

对于最低模态，仅有

最后，将所有角向频率的Galerkin系统排列在一起，可以得到全局分块对角系统

不同角向频率之间完全解耦，并且每个径向系统的左端矩阵都是对角矩阵。三角截断下的总模态数为。因此，在右端项的Zernike系数已经给定的情况下，全部径向系统的组装、求解以及向标准Zernike系数的转换都只需要次运算和的存储空间。

由此可见，对于圆盘上的Poisson方程，Zernike谱方法对应的线性系统本身具有非常简单的稀疏结构。实际计算的主要瓶颈并不是线性系统的求解，而是函数值与Zernike系数之间的正变换和逆变换。

Logan–Shepp Ridge Polynomials

Logan–Shepp多项式是定义在圆盘上的另一组正交多项式。与Zernike多项式按照角向频率与径向次数进行编号不同，Logan–Shepp多项式由总次数和方向指标编号，定义为

其中是第二类Chebyshev多项式。令以及，则。因此，每一个基函数只依赖于沿某个方向的投影，这也是它们被称为脊多项式（ridge polynomials）的原因。

Logan–Shepp多项式在单位圆盘上正交，并且次数不超过的三角截断为

截断空间的维数同样为。

Logan–Shepp基与Zernike基虽然具有完全不同的形式，但它们张成相同的有限维空间：两者的阶三角截断都等于笛卡尔坐标下总次数不超过的二元多项式空间。因此，在精确算术下，对相同函数使用相同截断次数时，两种展开给出的最佳多项式逼近完全相同，区别仅在于所采用的基底和计算方式。

对于固定的总次数，Logan–Shepp多项式与所有满足且的Zernike多项式之间存在离散Fourier变换关系。例如，当时，

奇数次数也具有类似的离散Fourier关系。这说明同一总次数下的Logan–Shepp基可以看成Zernike基经过一个离散角向变换后的结果。

Poisson方程的Dual Reciprocity方法

Logan–Shepp基对于Poisson方程的主要优势并不在于形成稀疏微分矩阵，而在于可以方便地使用Dual Reciprocity Method。由于只依赖于单一投影变量，寻找满足的特解可以转化为一维积分问题。

第二类Chebyshev多项式满足。因此，可以选择使得。具体地，

由于对应的方向向量具有单位长度，有。因此，

若右端项已经展开为，则一个特解可以直接写为

为了满足边界条件，再添加一个调和修正项

系数和就是边界函数的Fourier系数，可以通过FFT计算。最终解为。

该方法避免了对Poisson算子构造大型线性系统：计算过程被分解为右端项的Logan–Shepp展开、逐项生成特解以及一次边界Fourier修正。然而，Logan–Shepp基不是张量积基，无法像Zernike基那样进行径向与角向的分步求和。按照文中的直接二维求积实现，计算全部展开系数的代价可达到。因此，它适合右端项展开次数不太高，或者可以复用展开与特解信息的情形。

Chebyshev–Fourier Series

另一类自然选择是在角向使用Fourier级数，在径向使用Chebyshev多项式。文献讨论了以下四种径向基：

1. 在上使用完整的；
2. 在上使用平移Chebyshev多项式；
3. 使用二次变量变换后的；
4. 对角向频率使用与奇偶性相同的Chebyshev多项式。

由恒等式可知，第三种方法与偶数阶Chebyshev多项式等价。另一方面，圆盘函数的第个Fourier系数关于具有与相同的奇偶性，因此使用完整的会浪费约一半的自由度。文献最终推荐使用奇偶限制的Chebyshev–Fourier展开。

使用复Fourier形式，可以将这一展开简写为

这里偶数角频率只与偶数阶Chebyshev多项式相乘，奇数角频率只与奇数阶Chebyshev多项式相乘。该选择保证径向奇偶性正确，但并没有显式包含因子，因此一般的截断展开并不会自动满足原点处的全部解析性条件。特别是对于较大的，低阶Chebyshev项必须通过系数之间的抵消才能形成一个阶零点。

平移Chebyshev基会同时在和附近聚集节点。当函数的细小结构或复奇点靠近原点时，这种节点分布可能优于奇偶限制的基；但对于通常的圆盘问题，高梯度和近奇异结构更可能靠近外边界，因此奇偶限制的Chebyshev基一般更加有效。

Poisson方程的Chebyshev–Fourier离散

由于角向仍然使用Fourier级数，Poisson方程同样解耦为一系列径向方程

对每个，使用奇偶限制的Chebyshev基展开，并采用满足边界条件的Galerkin基，可以得到带状线性系统

Shen证明了相应的可以构造成七对角矩阵，即每一行至多包含七个非零元素。由于带宽不随截断次数增加，每个径向系统可以在复杂度下分解和求解，全部角频率的求解复杂度为。

Chebyshev–Fourier方法最重要的优势是函数值与系数之间的变换可以在两个方向上都使用FFT或离散余弦变换。当与具有相同数量级时，二维变换复杂度为，低于传统Zernike变换的。

其代价是自由度较多。若希望精确表示所有总次数不超过的笛卡尔多项式，Zernike三角截断只需要个基函数，而矩形Chebyshev–Fourier截断通常需要约个基函数。因此，对于一般非轴对称函数，Chebyshev–Fourier方法通常需要大约两倍于Zernike方法的自由度。

综合来看，在中小规模下，Zernike基通常以更少自由度达到同样精度，并且Poisson矩阵的带宽更小；当很大或需要在时间积分中反复进行正逆变换时，Chebyshev–Fourier方法的二维FFT优势会逐渐占据主导。

Cylindrical Robert Functions

圆柱Robert函数试图同时结合Chebyshev多项式的快速变换能力和圆盘函数在原点处的正确行为。按照文中的定义，其基函数为

并对施加相应的奇偶限制。因子显式保证了角向频率对应的径向系数在原点处具有阶零点，因此该基在结构上满足圆盘解析性的要求。

然而，这组基函数具有严重的数值病态性。当较大时，在除附近以外的整个区间上都极小，所有基函数都集中在宽度约为的边界层内。在这一狭窄区域中，低阶Chebyshev多项式变化非常缓慢，从而使、等函数几乎线性相关。

文献中的数值实验表明，固定每个取个径向函数时，插值矩阵的条件数会随着迅速增长，在较大角向频率下远远超过双精度浮点数所能处理的范围。因此，尽管Robert函数具有正确的原点行为并且理论上可以达到几何收敛，实际计算中却会因为基函数之间接近线性相关而丧失精度。

这说明将解析性因子直接乘在一个与无关的全局多项式基上并不一定是稳定的。Zernike多项式通过让Jacobi参数随变化，使径向振荡区域也随自然向边界移动，从而避免了Robert基的严重退化。Boyd和Yu因此认为圆柱Robert函数不应继续用于实际计算。

Fourier–Bessel Expansions

利用分离变量法求解圆盘上的Laplace或Poisson方程时，会自然得到Fourier–Bessel展开

其中表示Bessel函数的第个正零点。由于，所有基函数都自动满足处的齐次Dirichlet边界条件；同时，因此它们也具有正确的原点行为。

径向Bessel函数在权重下正交，

因此，第个Fourier系数的Bessel系数为

Fourier–Bessel基的最大优点是它由Laplace算子的特征函数组成。正确的特征关系为

因此，若右端项的第个Fourier分量展开为，则解系数满足完全对角的线性系统

即。从线性系统求解的角度看，这几乎是最理想的基底。

然而，Fourier–Bessel展开通常不是谱收敛的。原因是所有径向基函数及其反复作用Laplace算子后的结果都在处满足特殊的边界条件，而一般的光滑函数并不满足这些条件。对于固定，只有当

依次成立时，才可以在分部积分过程中不断消去边界项。

若前个兼容条件成立，即对成立，则相应Bessel系数的渐近衰减率为

每增加一个边界兼容条件，系数衰减阶数提高。只有当所有阶兼容条件都成立时，才可能得到指数衰减，而这对于一般解析函数几乎不会发生。

因此，Fourier–Bessel方法虽然具有对角的Poisson算子和自动满足的边界条件，但通常只能获得代数收敛。它最适合解析分离变量解、特征值问题或者右端项本身满足全部兼容条件的特殊问题，并不适合作为一般光滑圆盘函数的高精度变换基。

Square-to-Disk Conformal Mapping

另一种方法是通过共形映射将单位正方形变换到单位圆盘，然后在正方形上使用二维Chebyshev展开。文献使用余弦双纽线函数（cosine-lemniscate function）构造映射。若表示正方形上的复变量，表示圆盘上的复变量，则映射写为

其中是模数对应的第一类完全椭圆积分。

在正方形上选取张量积Chebyshev节点，再通过将其映射到圆盘。对于圆盘函数，定义拉回函数，并使用二维Chebyshev展开

由于正方形上的网格和基函数都是张量积结构，系数可以通过二维离散余弦变换在复杂度下计算。

共形映射不会破坏圆盘内部解析函数的谱收敛性。映射函数的复奇点位于正方形之外，因此拉回函数在正方形上仍然解析，Chebyshev系数保持几何衰减。不过，即使原函数是整函数，映射自身的有限距离奇点也会将原本可能出现的超几何收敛降低为普通的几何收敛。

该方法的主要问题是节点分布严重不均匀。正方形四个角附近的大量Chebyshev节点被映射到圆周上的四个特殊位置附近，形成四个非常密集的节点簇。除非解恰好在这些位置具有细小结构，否则这种聚集会造成明显的自由度浪费。

Boyd和Yu没有在文中详细推导这种方法对应的Poisson矩阵。由共形映射的标准性质，若和，则圆盘上的方程可以转换为

因此，使用二维Chebyshev配置法时，内部节点上的离散系统可以写成

并将边界节点对应的行替换为Dirichlet边界条件。这个矩阵具有Kronecker和张量积结构，但映射系数使其不再是简单的常系数可分离问题。

该方法适合已经具有正方形Chebyshev程序、希望以较小编程代价将其扩展到圆盘的情形。它能够保持谱精度，但通常不如直接使用适合圆盘几何的Zernike或Chebyshev–Fourier基高效。

Radial Basis Functions

径向基函数方法直接在圆盘上选取一组中心，并使用到这些中心的欧氏距离构造逼近。一般形式为

Boyd和Yu主要考察Gaussian RBF，

其中是形状参数，表示局部或全局节点尺度。

若在节点上进行插值，则系数满足

RBF方法不要求节点形成张量积网格，因此可以使用不规则节点，也可以在函数具有局部细小结构的位置自适应增加节点。这一灵活性是RBF方法相对于全局正交多项式方法的主要优势。

但是，RBF精度对节点分布和形状参数非常敏感。文献对个节点进行的实验表明，相比准均匀节点，将节点随机扰动通常会使误差增大约一个数量级，极端情况下可以恶化约三个数量级。在文中测试的三种节点尺度定义中，使用整个圆盘上的平均间距通常给出最好的结果。

对于Poisson方程，可以在内部节点配置微分方程，在边界节点配置边界条件。对二维Gaussian RBF，若，则

因此，Poisson配置系统可以写成

其中，而。

RBF矩阵通常是稠密的，并且不具有张量积分解。若总节点数为，直接分解需要次运算和存储；若用每个方向上的分辨率表示并取，则复杂度分别变为和。快速求和和迭代方法可以降低渐近复杂度，但实现复杂，而且常数通常明显大于FFT。

Gaussian RBF在形状参数趋于零的平坦极限下会逼近相应的多项式谱插值，因此经过良好调参后，其精度通常与Chebyshev–Fourier方法相当，而不是系统性地优于传统谱方法。与此同时，平坦极限中的插值矩阵会变得高度病态，需要RBF–QR等稳定算法才能可靠计算。

因此，对于规则圆盘区域，RBF通常比张量积谱方法昂贵；但对于非圆边界、不规则节点、局部自适应或难以构造全局坐标映射的问题，RBF提供了一种编程较简单且几何灵活的选择。

各种方法的比较

文献对上述方法的主要结论可以概括如下：

方法	一般光滑函数的收敛性	Poisson算子结构	主要变换代价	主要特点
Zernike	几何收敛	每个角频率对应低带宽矩阵	传统实现为	自由度少，适合中小规模
Logan–Shepp	与Zernike相同	可使用Dual Reciprocity避免直接离散微分算子	直接二维求积约为	Poisson特解容易构造，但非张量积
Chebyshev–Fourier	几何收敛	每个角频率对应七对角矩阵		可在两个方向使用FFT，但自由度约为Zernike的两倍
Cylindrical Robert	理论上几何收敛	可以进行谱离散	数值上无实际意义	严重病态，不建议使用
Fourier–Bessel	通常只有代数收敛	完全对角	径向变换不能直接使用FFT	Laplace特征基，适合分离变量和特殊兼容数据
Square-to-disk	几何收敛	正方形上的二维Chebyshev系统		易复用正方形程序，但四处节点过度聚集
RBF	调参良好时可与谱方法相当	稠密非结构矩阵	直接法为	几何灵活，适合不规则区域和自适应节点

Boyd和Yu的总体建议是：对于小到中等截断次数，Zernike基通常能够以较少自由度获得较高精度；对于非常大的截断次数或需要频繁正逆变换的时间依赖问题，奇偶限制的Chebyshev–Fourier方法更具优势。Logan–Shepp、Fourier–Bessel和共形映射方法分别适合特定问题，而圆柱Robert函数由于严重病态性应当避免使用。RBF的主要优势则体现在非圆区域和不规则节点上，而不是标准单位圆盘问题。