Computation on disk and sphere

REFERENCES

1. Computing with functions in spherical and polar geometries I. The sphere by Alex Townsend, Heather Wilber, and Grady B. Wright, SISC 2016.
2. Computing with functions in spherical and polar geometries II. The disk by Heather Wilber, Alex Townsend, and Grady B. Wright, SISC 2017.
3. Tensor calculus in polar coordinates using Jacobi polynomials by Geoffrey M. Vasil, Keaton J. Burns, Daniel Lecoanet, Sheehan Olver, Benjamin P. Brown, Jeffrey S. Oishi, JCP 2016.
4. Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations by Geoffrey M. Vasil, Daniel Lecoanet, Keaton J. Burns, Jeffrey S. Oishi, Benjamin P. Brown, JCP 2019.
5. Tensor calculus in spherical coordinates using Jacobi polynomials. Part-II: Implementation and examples by Daniel Lecoanet, Geoffrey M. Vasil, Keaton J. Burns, Benjamin P. Brown, Jeffrey S. Oishi, JCP 2019.

圆盘与球面是科学计算中两类重要的规则几何区域，广泛出现在流体力学、地球物理以及偏微分方程的数值求解中。与矩形区域上的张量积谱方法类似，我们希望利用 Fourier 多项式、Chebyshev 多项式或其他正交多项式构造高精度的函数逼近，并进一步实现高效的函数运算和 PDE 求解。然而，圆盘与球面的坐标结构使得这一问题并不完全等同于标准矩形上的谱逼近。例如，极坐标在圆心处、球坐标在两极处存在坐标奇异性，尽管函数本身可能非常光滑，其在坐标变量下的表示仍然需要满足特定的相容条件。因此，如何在保持函数正则性的同时，构造高效的谱变换与微分算子，是圆盘和球面谱方法设计中的核心问题。

针对这些困难，现有工作大致可以分为两条具有代表性的技术路线：

• 基于函数低秩结构的谱方法（DFS–BMC 路线）：Townsend、Wilber 和 Wright 的两篇工作分别针对球面与圆盘，通过适当的坐标延拓，将原区域上的函数转化为具有特定对称性的周期函数，并利用双重覆盖（double Fourier sphere/disk）及 BMC（block-mirror-centrosymmetric）结构保证延拓的一致性。在此基础上，作者使用连续 Gaussian elimination 自适应地构造函数的低秩分离表示，并结合 Fourier、Chebyshev 和 ultraspherical 谱方法，实现高效的函数逼近、微分、积分及部分 PDE 的快速求解。这条路线的主要特点是利用函数本身的低秩结构降低计算和存储成本，而不是首先构造适应几何的多元正交多项式基底。

• 基于 Jacobi 多项式的稀疏谱方法（Jacobi 路线）：Vasil、Lecoanet 等人的三篇工作则从另一角度出发，通过构造与极坐标或球坐标正则性相适应的 Jacobi 多项式基底，将坐标奇异点附近的光滑性要求直接包含在基函数中。其核心思想是利用 Jacobi 多项式的微分、递推和参数转换关系，使微分算子在不同参数的正交基底之间具有稀疏表示。对于球坐标中的向量和张量问题，作者进一步引入自旋加权球谐函数和适当的正交变换，以处理不同张量分量之间的耦合。这条路线更侧重于构造具有明确正则性和稀疏算子结构的谱离散方法，从而高效求解涉及标量、向量及张量场的 PDE。

这里必须先区分几何范围。文献[1]研究二维球面；文献[2]和[3]研究二维圆盘；文献[4]先建立上的张量微积分，再进入三维单位球；文献[5]主要实现和检验内的方程，也包含一个纯球面的 Rossby 波问题。因此，后三篇中的“spherical coordinates”不能统一理解成只在球面上求解方程。

下文逐篇沿用原文的局部符号。文献[1]用表示经度和余纬；文献[2]用表示角度和半径；文献[3]用；文献[4]、[5]用表示半径、余纬和经度。相同字母在不同论文中未必含义相同，尤其是低秩阶数、Fourier 模态、Jacobi 参数与球谐次数。

需要注意的是，前两篇文章从如何方便地计算一个函数出发：先利用几何对称性将函数延拓到规则矩形，再用保结构的低秩逼近把二维计算尽可能化为一维计算。它们的代表性成果是 Spherefun 和 Diskfun，以及与这些函数表示相容的快速 Poisson 求解器；后三篇文章从如何准确而稀疏地表示微分算子出发：先分析光滑标量、向量和张量在坐标奇点附近必须满足的结构，再构造直接满足这些结构的 Jacobi 基底，并允许微分前后的基底参数发生变化。这样，原本包含、以及分量耦合的表达式，可以转化为系数空间中简单的升降阶运算。因此，两条路线并非针对完全相同的计算问题。总体而言，DFS–BMC 路线更强调自适应函数表示及其快速运算，Jacobi 路线则更强调几何正则性与微分算子的稀疏性。下面分别介绍这五篇工作的主要思想、数学工具及具体计算流程，最后比较两条路线在逼近效率、数值稳定性和 PDE 求解方面的优势与局限。

DFS–BMC 路线

Computing with Functions in Spherical and Polar Geometries I. The Sphere

1. 从球面函数到双周期函数

这篇文章希望把 Chebfun 中“直接对函数进行数值计算”的思想推广到球面。目标不仅是插值，还包括求值、积分、微分、向量分析和 PDE 求解。其基本工具是 double Fourier sphere，简称 DFS，以及连续版本的 Gaussian elimination，简称 GE。[1, §1–2]

Step 1 用经度与余纬描述球面。 原文采用

其中。同一个符号也用于表示拉回参数域的函数：

经度方向天然周期，但余纬方向只有。如果直接在这个矩形上使用普通张量积插值，既没有利用越过极点后的连续几何关系，也容易在两极附近采集大量重复信息。

Step 2 通过滑移反射把余纬区间加倍。 将角度按周期理解，可以把原文式(2.2)等价地写为

这样得到的在两个方向都以为周期，并满足

这个关系有直接的几何意义：越过北极时，余纬反号与经度移动描述的是同一个球面位置。DFS 将这种关系变成可以由 Fourier 表示保存的代数对称性。

在上，令、。原文将此延拓画成分块形式

并称其为 block-mirror-centrosymmetric，简称 BMC。这里表示第二坐标方向的反射。球面还要求和处的函数值分别与经度无关；满足这些极点常值条件的 BMC 函数称为 BMC-I。南极两条边和还必须取同一个值。[1, Definition 2.1–2.2]

DFS周期延拓：左图为未延拓之前的参数区域，黑色加粗线段上的所有坐标点对应球面的同一南极点，灰色加粗线段上的坐标点对应球面上的同一北极点；右图为DFS延拓之后得到的参数区域，其中相同颜色的有向线段上的点的值按照箭头方向从头到尾逐一对应相等。具体来讲，左下块和左上块经过上下颠倒再平移可以分别直接得到右上块和右下块（无需左右颠倒）。

从上面的图片可以看到，左侧图中的深灰色线段原本是矩阵的一条边界，直接在左侧进行标准矩阵区域上的Chebyshev逼近会在边界处大量布点产生浪费，经过DFS延拓之后，这条代表北极点的深灰色线段变到了一个更大的区域的内部，这样一来消除了北极点处的大量重复节点，二来更好地保证了北极点附近的连续性和光滑性。

2. 为什么普通低秩分解需要修改

作者希望得到如下形式的低秩逼近：

这里和分别是一维纬向和经向函数，是尺度因子，是分离秩。只保存这些一维函数，比保存完整二维网格或系数矩阵更省空间，也方便复用一维算法。

普通连续 GE 每次选一个主元，减去一个行切片与列切片的乘积。但单个主元没有同时照顾几何上相配的点，因而一般会破坏 BMC 结构。作者把互相对应的点一起选入，形成一个主元矩阵。[1, §3.1]

[1, Figure 2]

对当前残差，令

一次更新为

其中是截断伪逆：只对大于阈值的奇异值取倒数。由于的奇异值就是与，这一步只涉及一个极小矩阵，并不需要对完整二维采样矩阵做 SVD。原文用控制两个对称通道的耦合，实验取。

这样设计同时解决了两个问题。第一，及其截断伪逆与交换矩阵对易，所以更新保留 BMC 对称性。第二，当函数在某个对称通道上很弱时，不会为了强行求逆而除以一个很小的主元；更新可以自动选择秩一或秩二。理论分析将该过程分解为两个有关联的标量 GE 过程：一个通道在余纬上为偶函数、在经度上为周期；另一个在余纬上为奇函数、在经度上为反周期。[1, §3.2–3.3]

极点的非零常值需要预先单独处理。原文式(3.7)的意图是先抽取一个与经度无关的分离项，使其余列在两极为零。严格保结构时，这个首项的纬向部分还应为偶函数；例如，可取某一列与其反射列的平均

它保留原函数的两个极点值，减去后仍为 BMC-I。这里写出的偶部投影，是把原文首项处理的保对称要求显式补齐；任意单列本身未必是偶函数，不能仅凭“减去一列”就推断 BMC 得到保持。

3. 自适应构造、理论保证与实际成本

实际算法先在较粗网格上寻找较大的块主元，估计所需秩与主元位置；随后分别加密行、列切片，以一维 Fourier 系数的衰减判断分辨率是否足够。利用 BMC 对称性，不必真的计算整个加倍区域。最终得到的函数表示既在两个方向上自适应，也对分离秩自适应。[1, §3.1]

设需要个 Fourier 系数，需要个 Fourier 系数。原文给出的构造成本为

保存低秩表示需要个数。这些优势都以明显小于为前提；“函数光滑”并不等于“在任意分辨率和坐标方向下秩都保持常数”。

理论上，文献[1]证明了理想化算法对有限秩 BMC 函数的精确恢复，并在足够大的复平面 stadium 解析区域等假设下证明残差一致趋于零。这里的 stadium 是围绕实区间的一条复平面邻域，两端由半圆封闭。解析区域的大小需要压过 GE 更新可能产生的增长。这篇球面文章把几何收敛率列为进一步的预期；下一篇圆盘文章才给出相应的几何率定理。因此，不能将本篇结论写成“对任意光滑函数都已证明几何收敛”。此外，理论中的连续主元搜索与实现中的粗网格搜索也应区别开来。[1, Theorem 3.2、Theorem 3.5 及其后讨论]

4. 低秩表示如何支持函数计算

求值与系数变换。 一次求值可以先分别计算和，再做项求和。需要二维 Fourier 系数时，先对每条一维切片做 FFT，再用其系数向量的外积组合；若确实要输出整个矩阵，就仍有写出这些元素的成本。原文分别讨论了保留因子和形成完整矩阵两种接口。[1, §4, pp. C416–C417]

积分。 球面面积元为，因此

积分化为一维加权积分与周期积分。如果，第二个积分为。这个来自原文采用的 Fourier 归一化，不能省去。[1, §4, pp. C417–C418]

微分。 原文主要通过嵌入后的笛卡尔切向导数组织球面向量分析。例如

这里沿用原文的切向导数含义，第三式的负号由确定。对一维 Fourier 因子求导只需乘波数，乘、只涉及相邻系数。对于表面上出现的，作者在系数空间利用表示“乘”的带状矩阵求解除法，避免在极点直接逐点除以零。原文在这一步采用偶数阶截断，使相应的有限乘法矩阵可逆。[1, §4, pp. C418–C419]

进一步将这些导数组合，就得到表面梯度、散度、旋度和涡度等操作。向量场由笛卡尔分量组成，各分量仍可用 Spherefun 表示。因此，第一条路线并非只能处理标量；它采用的是不同于自旋张量基的组织方式。[1, §4.1, pp. C419–C420]

5. 球面 Poisson 方程为什么可以快速求解

考虑

第一个积分条件保证方程可解，第二个固定常数零空间。球面没有外边界，因此这里的约束不能与圆盘边界上的 Dirichlet 条件混为一谈。[1, §5]

将 Laplace–Beltrami 方程乘以，得到

设是的双 Fourier 系数矩阵，为一阶 Fourier 微分矩阵，表示乘。原文式(5.5)为

其中表示的系数，而不是未经加权的的系数。令

由于，每个经向 Fourier 模态满足

这里的负号直接由二阶 Fourier 微分确定。时是固定带宽系统；时加入均值约束以消除常数零空间。于是，在右端系数已给定后，所有模态的求解总成本为。[1, 式(5.5)–(5.7)]

这个“最优复杂度”有明确口径：它针对求出个系数的线性求解阶段。若从完整网格构造右端，通常还有的二维 FFT；若从秩因子形成完整右端，外积装配一般为，把它记成隐含了有界。完整求解器输出通常也不是自动保持原输入秩的低秩函数。

原文报告在当时的标准笔记本上，约一分钟可以完成一亿自由度的 Poisson 求解。它说明这种带状分模态算法具有很低的实际常数，但不构成对一般变系数、非线性或张量 PDE 的统一复杂度结论。[1, §5、Fig. 7]

6. 本篇的主要贡献

这篇文章的贡献是把经典 DFS 延拓、自适应一维函数表示和保结构连续 GE 结合成一套可计算的球面函数系统。特别值得保留的思想有两个：几何对称性必须在压缩过程中保存；一个合适的矩形表示既可以支持低秩函数计算，也可以支持固定带宽的 PDE 离散。

同时，BMC-I 主要规定周期、对称和极点常值，不能单凭这些条件就推出任意近似函数在球面上的所有阶导数都光滑。对于光滑目标，文中给出了有效的自适应逼近和微分实验；这与“整个 BMC-I 函数空间自动满足全部极点正则性”是两个不同命题。这个区别会在两条路线的比较中进一步展开。

Computing with Functions in Spherical and Polar Geometries II. The Disk

1. 从双 Fourier 球面方法到圆盘的负半径延拓

第二篇沿用前篇“几何延拓与保结构低秩逼近”的框架，但圆盘的半径方向不应周期化。原文采用

利用与对应同一点的事实，将函数延拓到，使其满足

这是 DFS 在圆盘上的对应形式。角向仍用 Fourier，径向则在对称区间上使用 Chebyshev 多项式。与在上直接使用单边 Chebyshev 网格相比，负半径延拓把圆心放到区间内部，减轻节点在圆心的集中；保结构低秩采样又进一步减少完整张量网格中的冗余。[2, §2.1–2.4]

在 BMC 对称性之外，圆盘要求与无关。原文把这一条件称为 BMC-II。若写成角向 Fourier 展开

则对称性和原点单值分别给出

所以偶角频对应偶径向函数，奇角频对应奇径向函数。但这些条件尚未要求在原点具有阶零点；后者是物理解析性提出的更强要求。[2, Definition 2.1、§3.3]

2. 保结构低秩逼近与新增的收敛分析

圆盘的 GE 主元与前篇相同，只需将替换为。每次使用一个中心对称的主元，通过截断伪逆保留一个或两个奇偶通道。原点处的非零常值先由角向常数项承载，其余非恒定角向项的径向因子在原点为零。[2, §3.1–3.3]

原文式(3.10)仍把首项写为任意径向切片。与球面情形相同，严格保 BMC 时应取该切片的偶部；如果只需归零圆心，直接减去原点常数也是合法的保结构选择。这是对几何约束的明确解释，而不是断言任意切片都自动具有正确奇偶性。

最终表示为

其中

这里是径向系数数，是角向系数数，不能与第三篇中的 Fourier 阶混用。若分别收集这些一维系数，，则完整 Chebyshev–Fourier 系数矩阵为

这里是普通转置，不是共轭转置，因为原分离表示是两个函数的乘积。构造费用仍为，样本量约为；取得低秩系数因子还需要相应的一维 FFT/DCT。若要显式形成，另需。[2, §3.2、§4.1–4.2]

本篇比球面篇进一步给出几何收敛率。原文 Theorem 3.4 考虑实值 BMC 函数：一个变量连续，另一个变量可在足够大的 stadium 内解析延拓并一致有界。按文中的角向解析假设，若

则理想化 GE 近似满足

取时。这里是解析区域的裕量参数，是 GE 步数，都与后面 Helmholtz 波数或 Fourier 模态的局部记号不同。证明的核心是控制两个奇偶通道之间的增长，使解析逼近带来的衰减占优。原文也明确承认该解析区域假设较强；实际算法对更一般函数的良好表现主要来自数值实验。[2, §3.4]

文中还将 GE 与连续 SVD 的最佳秩逼近进行比较：解析例子中两者显示几何衰减，有限光滑例子中显示代数衰减，GE 的误差曲线接近 SVD。这里的“接近最优”是这些实验的结论，且对应延拓矩形上的标准范数，不是对所有函数的统一近最优性定理。[2, §3.5、Fig. 5]

3. 函数运算与圆盘面积权下的低秩分解

求值与积分。 径向 Chebyshev 因子可用 Clenshaw 算法评价，角向 Fourier 因子可用 Horner 形式评价。面积积分则分解为

其中是面积元带来的权，不能用无权矩形积分代替。角向积分只取零频，径向积分可复用一维 Chebyshev 求积技术。[2, §4.1–4.3]

笛卡尔微分。 原文使用

径向导数用 Chebyshev 系数递推，角向导数和三角乘法用 Fourier 系数完成。对于除以的步骤，先构造表示“乘”的三对角矩阵，再在系数空间求解商函数，而非在圆心直接相除。原文在该除法步骤采用偶数阶截断，以保证有限乘法矩阵可逆。非恒定角向项在圆心归零，为这种多项式可除性提供条件；承载原点值的首项在角向求导后为零。[2, §4.4]

这种实现避免了显式的，但不应把“可以在系数空间完成一次除”理解成“已经强制了所有角频所需的全部消失阶”。对高阶物理微分，仍需区分这些不同层次的正则性。

物理面积权下的 SVD。 这是本篇相对于球面篇的一项扩展。为了计算真正的圆盘范数，作者使用

对已有的径向与角向低秩因子分别做相应内积下的 QR，再对小的核心矩阵做 SVD，便得到

径向 QR 的构造中用到了关于正交的 Bessel 参考函数。这里并不是把 Diskfun 的主要表示改成 Bessel 级数；主表示仍为 Chebyshev–Fourier 低秩因子。面积权 SVD 的最优截断也未必保留原点连续性，所以它用于范数等运算，而不能直接替代 BMC-II 构造。[2, §4.5]

向量场。 Diskfunv 以两个笛卡尔分量表示平面向量场，再通过上述导数形成梯度、散度和旋度。这样避免把原点处未定义的极坐标分量直接当成普通光滑标量。例如，常向量在极坐标中的分量依赖，但它作为笛卡尔向量场完全光滑。[2, §4.6]

4. 圆盘 Poisson 求解器的完整设计

本篇求解

正文详细推导的情形。非齐次边界可先作提升：选取光滑使，写，再求与齐次边界。这是对原文脚注中“归约为齐次问题”的标准解释，原文没有指定某一个具体的提升函数。[2, §5]

Step 1 消去显式分母。 在加倍区域上将方程乘：

Step 2 把外边界条件放进试探函数。 令

于是处自动为零。这里未知矩阵表示；若表示 Chebyshev 基下乘的矩阵，则真正解的系数是。

代入后，径向链式法则给出

前三项构成原文的径向算子。其中和来自对边界因子的求导，不能把它们漏掉。[2, 式(5.1)–(5.3)]

Step 3 选择稀疏的输出基底。 原文利用

将这一二阶操作表示到第二类 Chebyshev 多项式基中。再结合以及低次多项式乘法，得到固定带宽的径向矩阵。这一步已经使用了 ultraspherical 谱方法“让导数落在另一组基底”的思想，所以两条路线在算子离散层面并非完全分离。[2, 式(5.4)]

记为系数到系数的转换，为乘，为角向二阶微分。离散方程为

此处是本身的 Chebyshev–Fourier 系数，与球面篇中已含权的不同。[2, 式(5.5)]

Step 4 写成适合低秩求解的 Sylvester 方程。 原文定义

这里的逆是数学表达，不要求真的形成稠密。例如，对任意移位，

因此 ADI 所需的移位求解仍可通过固定带宽矩阵完成；角向为对角阵，相关操作更便宜。这个等价变换说明了为何能同时利用 Sylvester 方程的谱结构与原离散的稀疏性。[2, 式(5.6)]

5. Factored ADI 的作用、费用与限制

与前篇主要求出完整系数不同，本篇希望直接返回一个低秩解，使其可继续作为 Diskfun 对象计算。若右端具有秩的分解，factored ADI 就只更新径向与角向因子，不形成完整未知矩阵。每轮完成若干移位带状求解，经过轮得到

中间因子宽度最多约为，随后进行压缩。这一是 ADI 因子分解的原文约定，不应与前面的混淆。[2, §5]

ADI 的效果依赖两个算子谱的分离。原文特别指出不是正规矩阵：若，误差界还要包含。文中的迭代数估计为

其中由两个谱区间决定。约随二次增长、随略快于二次增长，是作者的数值观察，不能替换成已经证明的全参数渐近界。

计算成本需要保留两部分：

若只写第一项，就漏掉了可能占主导的压缩工作。文中结合迭代数的增长写出带有因子的总体费用；它的优势属于右端秩较低、ADI 迭代数较少且解可压缩的情形。

当这些条件不理想时，同一个广义 Sylvester 方程也可按 Fourier 模态直接解耦，每个模态解一个带状径向系统，求全部系数的成本为。若还要将结果重新变成低秩 Diskfun，则需计入网格恢复和 GE 压缩。原文软件会在 ADI 不合算时选用这种全秩求解路线。[2, §5, pp. C259–C260]

原文计时中，在且右端秩小于的部分实验里，ADI 相对比较的全秩流程最多约快五倍；在、、右端秩的实验中，报告不足五秒得到低秩解。后一例的自由度被人为增大以展示渐近行为，且计时对应当时的 MATLAB 与笔记本环境。因此，这些结果支持“小秩大规模问题具有明显优势”，不支持“所有圆盘 Poisson 问题都由 ADI 更快解决”。[2, Fig. 7 及脚注]

6. 本篇的主要贡献

圆盘篇并非仅把球面代码换一套坐标。它增加了 BMC-II 的圆心处理、Chebyshev–Fourier 混合表示、更明确的 GE 几何收敛分析、物理面积权下的低秩范数计算，以及直接输出低秩解的 Poisson–ADI 流程。整体上，它把函数构造、函数运算和 PDE 求解组织在同一种压缩表示中。

其限制也来自同一选择：效率依赖函数和解的分离秩；面积权下的最优压缩不自动保留圆心约束；一般变系数或非线性方程也不会自动继承本文 Poisson 算子的可分离结构。对这类问题，仍需另外分析秩增长、模态耦合与时间演化中的误差。

Jacobi 路线

Tensor calculus in polar coordinates using Jacobi polynomials

1. 核心问题是圆心的完整正则性

这篇文章研究圆盘上标量、向量和张量 PDE 的谱离散。它并不以获得一个低分离秩的函数对象为主要目标，而是希望直接从几何结构出发，使微分、径向乘法和边界操作都具有简单的系数表示。[3, §1]

Step 1 分离角变量。 沿用原文坐标，

角向微分成为乘。但在笛卡尔坐标下解析的标量，径向分量还必须具有

其中在原点附近解析。因而高角频不能在圆心附近任意取值，它必须相应高阶地消失。这比偶、奇性以及非零角频在原点取零更强。

Step 2 将这个结构直接写进基函数。 本节采用，其中表示非负有效径向阶；完整 Fourier 指标则属于。原文的径向公式主要按非负书写，完整 Fourier 级数的径向阶应相应理解为。定义

其中关于正交，原文的平方范数为

因此

注意，积分权中没有再乘，因为两个基函数已经各自带有。是关于的多项式次数，物理径向次数为；是基族与权函数参数，不是张量秩。对于，

正是适当归一化的 Zernike 径向族。[3, 式(4)、(7)–(8)、(15)–(16)]

这一步使任何有限截断都具有正确的圆心行为。但本文真正的新意还在下一步：不要求所有导数继续用原来的基底表示。

2. Jacobi 参数变化为何使微分只有一条非零带

对任意可微，链式法则给出

第一式中的恰好抵消的导数，输出反而具有更高的圆心消失阶；第二式则降低一个径向阶。这两种操作与 Jacobi 多项式的参数升降关系匹配。代入归一化后得到

第二式按理解；第一式在时为零。时两种物理微分组合相同，均按升到径向阶的关系处理，不能生成非法的。[3, 式(5)–(6)、(17)–(20)]

为了方便记录输入与输出基，原文使用 Dirac 记号。对一般数值工作者，只需把

看成一个由基函数组成的行向量，而是系数列。于是。矩阵在基函数恒等式中写在右侧，在作用系数时则从左侧相乘。

原文定义

从而

只有第一条上对角带，只有主对角带。稀疏性来自输出空间的改变；若把结果全部转回原来的基底，就不能继续用这两个单带矩阵概括微分成本。[3, 附录A, (A.1)–(A.4)]

3. 从标量推广到向量和任意阶张量

极坐标中的单位向量随旋转。对向量求导时，必须同时求导其系数和基向量；只给标量公式添加一个分量下标是不够的。

作者使用复旋转基

对一个标量 Fourier 分量，

这组基使极坐标连接系数对角化：

将对系数与对基向量的微分合并，就得到有效角阶。对秩为的张量，令

原文式(36)–(37)的展开和求导规则为

其中省略的参数由它的输入基决定。以上沿用原文在非负有效径向阶下的简写；若，径向参数须取绝对值，并按相应对称性折叠升降关系。完整整数 Fourier 指标下的正则因子是。这也解释了为什么同一 Fourier 模式的不同向量分量不能全部使用作为圆心因子。[3, §3]

将两次微分作张量收缩，可得标量 Laplace 算子

两条可能的升降路径满足

其中第二条路径需在允许的径向阶范围内解释。这个关系对应平面上的协变导数可交换。注意每一步的都随路径变化，并不是两张固定矩阵任意交换。[3, 式(39)–(41)]

Remark: 复基的正交性与张量收缩

原文的普通双线性收缩满足和，与复向量的 Hermitian 内积不是同一个运算。Laplace 算子前面的因子来自两种交叉收缩，不能用 Hermitian 正交关系把它们删掉。复单位变换有利于数值条件，并不意味着张量运算可以随意把转置替换成共轭转置。

4. 乘法、转换与 Heisenberg 代数

乘半径。 原文构造两个双带矩阵：

具有主对角与第一上对角，具有主对角与第一下对角。它们代表相同的物理乘法，但输出基不同，在对应空间之间互为转置。于是

表示乘，并且是对称三对角 Jacobi 矩阵。[3, §4, 式(42)–(49)]

转换基底。 同一个函数的参数可通过双带矩阵提升：

其非零元素为

不存在的项不予生成。转换改变的是同一函数的系数表示；例如，微分两次后的输出在基中，零阶项和右端就应使用转到该空间后再相加。[3, §5、(A.11)–(A.13)]

轴对称变系数。 若径向系数写成，其系数乘法矩阵就是。将近似为次数的多项式级数，

若所选多项式满足

原文用矩阵 Clenshaw 算法反向递推：

每次乘三对角至多使半带宽增加一。因此，次数的多项式系数给出总带宽至多的乘法矩阵。这种稀疏性以系数的逼近次数较低为前提；复杂系数会扩大带宽，非轴对称系数还会耦合不同，后者不在本文线性算子的主要研究范围内。[3, 式(50)–(56)]

从实现角度看，首先是无限乘法算子的多项式。有限矩阵计算时应保留足够的末端填充再截取目标块，避免边界截断改变需要的矩阵元素。这是上述算子恒等式的实现要求，不能由“使用 Clenshaw”自动省略。

这些算子进一步满足

并配有相应的交换关系。原文将它们组织为二维 Heisenberg Lie 代数；严格解释时，可将所有空间取直和。实际用途是给微分、乘法和换基规定一致的代数规则，使复杂 PDE 可以由少量基本算子自动组装。[3, §5.1]

5. 外边界条件如何进入离散系统

圆心正则性已由基底处理，不应再把当作独立的物理边界。真正需要施加的边界条件位于。原文给出

以及

因此 Dirichlet、Neumann 或 Robin 条件都可以组装成系数上的已知行向量。文中多个实际例子采用行替换：用边界迹替换内部方程最后一行；在广义特征值问题中，将质量矩阵对应行置零。这样内部矩阵带状，但边界行一般稠密，完整系统应称为 almost-banded。[3, §5.2、§6 Example 1]

作者还讨论了利用换基来局部化边界条件，其基础是正确的转置恒等式

所以对齐次 Dirichlet 问题，可以在乘以一组正则 Jacobi 基的空间中寻找解，系数变换只增加一条带。对于一般边界值，一个完整的边界适配多项式基可明确写为

其中首项在取，其余项取，且多项式次数依次为。这组带列移位的写法是由原文基础恒等式得到的明确补充，不能直接与原文式(67)的矩阵等同。

Remark: 边界转换公式的两点校读

第一，附录(A.15)将的第一下对角系数印成正号，但由正文的转换与乘关系，该项应为负。例如，

若将负号换为正号，右侧会成为。因此上文采用正文的正确转置恒等式。

第二，原文式(67)写，没有列移位。它形式上得到，但不能据此断言在所有光滑函数之间构成完整可逆换基。考虑：若表示光滑函数，则边界式要求其首系数为零，而除以后又要求恰为第零个 Jacobi 基函数，首系数应为一，产生矛盾。上文的通过在齐次部分采用索引补齐了这一问题。这一校核并不影响本文数值例子所用的显式边界迹与行替换方法。

对于 Neumann 条件，常数模态的边界导数为零，所以任何试图通过其边界导数系数取逆的公式都须排除该模态，并用均值或其他规范条件固定常数；方程本身也需满足相应的可解性条件。[3, 式(72)–(75)]

6. 从网格值到 PDE 解的实际流程

Step 1 角向 Fourier 变换与分量变换。 对共同径向网格上的函数做角向 FFT；对向量和张量，把坐标分量变换到复旋转基，并为每个分量使用正确的有效径向阶。

Step 2 在共同径向网格上投影。 原文取个 Gauss–Legendre 节点，令

若标准 Legendre 权重为，则径向权重，因为。合成与分析分别为

采用

并仅保留的模式后，被投影的有限展开在中的乘积次数不超过，因此上述求积对该有限多项式空间精确。实际通常先求系数，再由提升到其他。这里的“精确”是指被充分分辨的有限展开，不是说任意连续函数只用有限节点就可精确投影。[3, §5.3]

若，，则有限空间满足。分析使用加权转置，合成使用；未加权采样矩阵本身不必是方阵或正交矩阵。

Step 3 统一输入输出基后组装方程。 以

为例，对每个定义

这里是为说明流程引入的缩写。内部方程为

与都把输入表示到基中，随后加上。有主对角和两条上对角带；原文称它为 tridiagonal，准确的中文解释应是“三条非零带”，避免误认为通常的对称三对角结构。[3, 式(84)、(86)、(110)]

Step 4 求解与还原。 对轴对称线性系数，各独立，可使用带状或 almost-banded 线性代数，再作径向矩阵评价与逆 FFT。Example 4 使用 ApproxFun.jl 的 adaptive QR，按每个模态的需求决定径向长度。[3, §6]

7. 四个数值例子分别验证了什么

Bessel 特征值问题。 原文将与写成广义特征值问题，并与 Bessel 零点比较。在的例子中，前约个充分分辨的模态达到接近浮点精度；第个特征函数保持的圆心行为，图示误差约。高频欠分辨部分仍有大误差，不能概括成全部离散谱都精确。[3, Example 1、Figs. 2–4]

无粘惯性波。 对旋转圆柱的速度–压力系统，作者使用径向基表示速度，将梯度、散度和旋转作用组装为少量块，并施加壁面无穿透条件。算例检验了真正的向量分量耦合；所得频率仍位于物理区间，但部分欠分辨模态的误差显著。[3, Example 2]

Hagen–Poiseuille 管流稳定性。 背景轴向流为。流速乘法通过的组合表示，剪切项通过乘表示，黏性项通过两次表示，压力和不可压缩约束保留在原始变量块系统中。文中在不同 Reynolds 数下与已有特征值基准比较；这是轴对称背景的线性稳定性计算，不是非线性湍流模拟。[3, Example 3、Table 3]

强迫 Helmholtz 与 screened Poisson。 右端和边界数据包含完整角向依赖，但常系数算子仍可逐求解。随着波数增大，使用 adaptive QR 增加分辨率；原文报告最大测试超过万未知量、求解少于秒。该计时属于特定系数求解流程与当时环境，不应与其他论文不同问题的时间直接横向比较。[3, Example 4、Figs. 7–8]

8. 本篇的贡献与尚未解决的瓶颈

本篇的主要贡献是建立一套由标记的基底演算：复旋转基处理连接项，Jacobi 参数升降处理微分，转换矩阵统一不同导数阶的输出空间，径向乘法和边界迹支持完整 PDE 装配。它使向量和张量方程可以在原始变量下保持稀疏，而不必先改写成更高阶的纯标量方程。

不过，正反径向变换仍是稠密矩阵乘法。对一个，费用为；若径向和角向分辨率均与同阶，总径向变换通常为，另加角向 FFT。原文的线性求解费用不能转述为包含网格变换的整个圆盘算法复杂度。[3, §5.3、Example 4]

原文一维计时中，矩阵乘变换与 FCT 在阶附近相当，阶时前者约为后者的倍，阶时约为倍。这说明中小规模矩阵乘法可能有竞争力，同时也展示了渐近瓶颈。[3, Table 1]

最后，2016 年本文实际给出的四组实验均为线性、时间无关问题或特征值问题。文末讨论了时间演化的可能性，但没有给出非线性去混叠、具体时间推进器或长期模拟验证；这些内容应放到2019年的实现篇中介绍。

Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations

1. 从圆盘推广到球内，新增的困难是什么

这篇文章的最终目标是三维单位球内的张量 PDE。原文采用

体积元为。相比圆盘，既有南北极的角向退化，也有球心的径向退化。作者先在上建立角向张量微积分，再寻找能够将各分量的球心行为分开的变换，最后才构造径向 Jacobi 基。[4, §1]

标量球谐分量满足，但向量的三个分量不是简单地各带一个。这正是从二维圆盘到三维球内最关键的新问题：自旋基可以整理两极的行为，却还不能单独整理球心的行为。

2. 球面上的自旋基与自旋加权球谐函数

Step 1 对角化球面连接。 与圆盘类似，作者定义

对标量，、。在球面切空间内，

若是若干基向量的张量积，只需以总自旋替换单个。[4, §2.1, 式(19)–(25)]

Step 2 每一种自旋使用相应的球谐基。 原文将自旋加权球谐记为，其精确 Jacobi 表示为

这里

并采用原文的相位与归一化

时恢复普通球谐。因子控制北极行为，因子控制南极行为，两个指数同时依赖经度波数和自旋。因此，一个自旋非零分量如果仍按普通标量球谐处理，会丢掉本来应与基向量旋转相匹配的极点结构。[4, 式(26)–(28)]

Step 3 把协变微分变成自旋升降。 原文给出

其中

将它与基向量的导数合并，

完整梯度再乘并对求和，张量秩随之增加一。这说明“球面微分像 Fourier 微分”的确切含义：固定时只乘一个明确系数，但输出进入自旋改变后的另一组基，而不是始终停留在同一基中。[4, 式(29)–(34)]

Step 4 区分不同的 Laplace 算子。 两次内禀协变微分收缩后有

这是原文所称 connection 或 rough Laplacian 的系数。只有时直接等于。三维向量或张量 Laplacian 限制到球面时，还可能出现曲率相关的零阶差异；不同物理方程采用哪种算子，需要根据方程本身判断。[4, §2.3]

3. 角向乘法也有稀疏结构

除了微分，PDE 还需要乘、等几何系数。固定，原文给出

其中

或不存在的模式须按相应极限与零系数处理，不能直接计算形式上的。乘保持自旋、只耦合相邻球谐次数；乘则可以通过的三项关系表示。原文把这些矩阵记为，与后面的径向转换不是同一个矩阵。[4, §2.4, 式(40)–(49)]

例如，球面 Coriolis 项中的只引入少量耦合。一般变系数仍需按其展开长度组合乘法矩阵；简单系数的三项关系并不意味着任意变系数都具有固定带宽。

4. 为什么球内还需要一个依赖于 的正交变换

标量径向正则性空间在原文中记为

这里的沿用原文表示规定的局部结构，不宜理解成仅保留一个首项而允许任意奇异余项。它满足

在固定下，向量可以分解为

但这三个正则性分量并不等于三个自旋分量。作者通过分析向量 Laplacian 的角向耦合，发现对应特征值为、和，由此得到分离它们的正交变换。[4, §3.2]

令

对于，按顺序，原文式(75)为

将正则性系数转为自旋系数，则沿相反方向作用。由于它依赖，必须先做球谐变换，才能完成这一分量旋转；不能在物理网格上用一个固定旋转矩阵提前替代。

对于更高阶张量，问题还有一个新层次：不同分量可能具有相同的正则性，单纯对角化 Laplacian 不能唯一决定简并子空间内的基。例如二阶张量在充分大的下有九个分量，分属

作者因此利用梯度的递推结构构造一般的，通过表示矩阵的 Kronecker 和记录各个张量指标的作用。这使变换不但正交，而且与微分递推一致，而不是随意选择每个重复特征值的特征向量。[4, §3.3, 式(87)–(100)]

小需要删去不存在的自由度。时，相应自旋球谐不存在；例如的向量只剩对应的径向分量，不能保留形式上的。原文专门讨论了这类例外。[4, 式(101)–(103)]

5. 径向张量微分归结为两个算子

作者定义

它们的组合给出标量径向 Laplacian：

在相应通道存在时，把映到。这里将整个组合当作一个微分操作，而不是分别计算可能很大的导数项与项后再相消。[4, 式(51)–(55)]

设标量径向系数为，向量在正则性基中的系数为。经过变换后，梯度与散度变为

若，则

这些公式的形式很简单，但其中已是正则性分量，不能当作原来的坐标分量。[4, 式(76)–(86)]

更一般地，若张量系数写作，完整展开为

原文的任意阶梯度规则可写为

的输出为零，无需引入。这个公式说明高阶张量仍使用同一组径向算子，不必每增加一阶就重新推导所有球坐标连接项。[4, 式(104)–(116)]

6. 广义球 Zernike 基为什么需要第二参数

取，对有

第二式中的与 Jacobi 升降关系匹配，所以径向 Jacobi 的第二参数应取。这一步既反映三维体积元，也反映三维微分算子的具体系数；不能只把圆盘公式中的替换成。[4, 式(117)–(123)]

原文定义

广义二项式按 Gamma 函数理解。基底满足

对应物理径向测度，此时归一化因子简化为。仍是关于的次数，物理径向次数为。[4, 式(124)–(125)]

微分关系为

第一式降低，第二式保留，两式都将提升一。降正则性公式只在对应物理通道存在时使用；小情形应结合前面的分量删减处理。[4, 式(126)–(128)]

原文同时用表示分量旋转，用表示径向基函数，二者作用在不同指标上。下文的则是收集所有径向基函数的行向量。

7. 转换、乘法与系数方程

原文的径向矩阵恒等式为

和都是双带矩阵；系数微分矩阵为对角，只有第一条上对角带。应区分函数上的与系数上的。[4, 式(134)–(143)]

例如，转换基函数的两项关系为

其中

这个负号决定了转换转置对应乘的正确关系。

考虑原文示范的标量问题

固定并令后，Laplacian 落入基，零阶项需提升到同一检验基，于是内部方程为

它清楚展示了稀疏谱离散的原则：不同微分阶的项在共同输出空间汇合，而不是把所有导数拉回原试探基。[4, §5.1, 式(151)–(156)]

8. 修正与 Galerkin 边界基的区别

方法一：用一个预定形状的残差承载截断误差。 原文在内部方程加入，其中的形状预先指定、未知，并另外施加边界迹。的物理径向次数为；本处以编号，共有个径向系数。在有限系数空间中，得到

是残差多项式在基中的系数。若取，则，为对应最高保留的第阶 Jacobi 模态的单位向量，这与替换最后一条内部方程等价。若取，则

只有三个可能非零的位置。两个系统都很稀疏，但允许承载残差的多项式不同，因此有限分辨率解和误差可以不同。[4, 式(157)–(162)、Appendix B]

方法二：直接让试探空间满足边界条件。 利用前面的转置关系，令

则

因此自动成立。把内部矩阵右乘，就得到对新系数的带状方程；对文中这个二阶标量问题，非零带只向下增加一条，不再需要稠密边界行。这是 Galerkin basis recombination。有限截断时须保留乘法产生的最高阶并匹配方程数，不能先截断乘法再假定边界恒等式仍严格成立。[4, 式(163)–(165)]

作者强调的一个重要观点是：算子的稀疏实现与截断残差的选择可以分别设计。 仅仅更换表示同一有限残差问题的基底，不会改变精确算术下的解；改变多项式则可能改变离散问题本身。下一篇会用实际算例说明这一差异。

9. 本篇贡献及其计算边界

本篇方法的核心结构有三层：自旋球谐统一球面角向张量微分；依赖的正交变换分离球心正则性；广义球 Zernike–Jacobi 基使径向微分和简单乘法成为单带或双带运算。角向自旋框架承接已有工作；本文进一步给出适合计算的角向乘法关系，并重点发展球内正则性变换与径向算子体系。对于任意阶张量，还通过递推避免重新手工推导全部坐标公式。[4, §2、§3、§6]

这些是精确的几何与算子结构，并不等于已经获得了快速自旋球谐或快速径向 Jacobi 变换。原文§5.5讨论了快速连接算法的可能性，也指出逐次改变参数在与分辨率同阶时仍可能昂贵。因此，不能把这篇文章关于稀疏算子的结论改写成“已实现接近 FFT 成本的完整球内变换”。

同样，正则性适配避免的是一类由错误基底带来的消去和奇点问题，并不自动给出所有 PDE 的条件数、所有参数范围的舍入误差界或所有离散特征值的稳定性保证。完整计算效果需要结合下一篇的变换、边界与时间离散实验判断。

Tensor calculus in spherical coordinates using Jacobi polynomials. Part-II: Implementation and examples

1. 从算子公式到完整求解器

Part-I回答“怎样表示和微分”，Part-II进一步回答“怎样变换、装配、施加边界、推进时间并验证结果”。作者利用 Dedalus 的 Fourier 变换和 MPI 数据转置等基础设施实现球内算法，依次检验特征值问题、边值问题以及非线性流体与磁流体问题。论文结论仍把从方程字符串自动装配全部矩阵和非线性项列作后续工作，因此这里讨论的是当时论文中的具体实现。[5, §1、§7]

原文把一般初值问题写为

是由不同阶张量组成的状态向量，是线性算子，可以包含非线性运算。这里表示边界函数；后面对流问题中的表示 Ekman 数，径向基转换始终记为。[5, 式(1)–(3)]

单个张量场的统一展开为

其中是球坐标分量的多重指标，是自旋多重指标，是正则性多重指标，上横线表示各指标之和。整理坐标分量与自旋分量，整理自旋与球心正则性，才是径向基函数。原文相近的三个记号在实现中对应三个不同步骤，不能合并成一次普通球谐变换。[5, 式(4)–(5)、Table 1]

截断首先受到几何结构约束。 对实场，程序利用共轭关系只保存；相应自旋的球谐次数满足

径向多项式总次数为。按照原文§2的设计，若计算需要的最高张量阶数为，则采用

只保留具有合法整数的模式。使用共同的，可以使同一下各正则性分量的径向系数等长，简化方阵装配。这种三角截断与“每个都保存相同数量径向多项式”的矩形截断不同；比较其他代码的分辨率时，应比较实际自由度与可表示的次数，不能仅比较这一名称。[5, §2、§6.1.2]

Remark：原文截断范围的一处不一致。

本节采用§2正文及 Table 1的总次数限制。§3在式(7)之后另写，它与前述范围并不等价；不能将二者当作同一公式使用。类似地，§3对的一处表达缺少总自旋的绝对值，这里采用§2明确给出的。这些说明只涉及原文两处计数表达，不能据此推断全部数值代码采用了哪一种未经核对的截断。

2. 网格变换的六个步骤

Step 1 建立公共求积网格。 经度使用均匀网格，余纬使用 Legendre–Gauss 求积，径向则在上使用参数的 Jacobi–Gauss 求积。所有的径向基在同一组物理点上评价，这使不同模式能够还原到共同网格上计算非线性乘积。[5, §3]

径向权重与物理测度的关系是

因此，若是权对应的 Jacobi 权重，则物理径向权重为。这里选来自共同的体积测度，不能误写成对每个分别使用的不同节点。

Step 2 预计算基函数与变换矩阵。 作者用 Golub–Welsch 算法给出求积点和权重的初值，再作 Newton 修正；基函数由端点包络乘 Jacobi 三项递推得到。高时，可能先下溢到零，后续递推便无法恢复本来应当具有中等大小的基函数值。实现因此在网格、权重和矩阵的构造阶段使用128位精度，之后转成64位保存并用于反复变换。论文报告这种处理在约以内的多项式次数上有效；这属于具体实现经验，并非任意规模下的稳定性定理。[5, §3]

Step 3 做经度 FFT 和自旋角向变换。 从球坐标分量出发，先对做实到复 FFT，再通过转为自旋分量。对于每个、每种自旋和每个径向点，用预存的加权球谐矩阵求系数。MPI 转置在需要时使当前变换方向的数据完整地保存在本地进程中；它改变数据分布，不改变数学表示。[5, §3, 式(6)]

Step 4 作用依赖于的正则性变换。 角向变换完成后，再用从自旋系数进入正则性系数：

这个次序不可颠倒，因为在球谐变换前，物理网格数据还没有一个确定的可供选择。[4, 式(107)；5, §3]

Step 5 做径向矩阵变换。 状态变量先展开在基中：

同一和正则性分量的径向矩阵可以在不同之间复用。逆变换依次执行径向评价、、自旋球谐评价、坐标分量恢复和逆 FFT。若直接投影到一般基，内积还须包含；实际状态系数从出发，各项通常通过提升到需要的输出基。合成矩阵与分析矩阵的关系是求积意义下的加权转置，不能去掉权重后将二者称为普通矩阵互逆。[5, §2–§3, 式(7)；4, 式(125)、(134)]

Step 6 在系数空间微分，在物理空间相乘。 例如计算，先以稀疏算子求梯度，再把速度与梯度评价到共同网格，逐点进行张量缩并，最后分析回谱空间并将各项转换到方程的共同输出基。这样能够复用 Part-I 的任意阶张量微积分，而无需显式装配一个巨大的非线性模式卷积张量。[5, §2、Appendices B–D]

原文同时报告了去混叠和不去混叠的实验。其网格规模为

设置	经度点数	余纬点数	径向点数
不去混叠
去混叠

需要时将点数取为合适整数。已充分分辨的那些算例没有显示去混叠带来系统性的精度差异；这不能推出一般非线性问题无需去混叠。[5, §3]

3. 线性装配、时间推进及真实计算成本

在本文主要三维算例中，隐式线性算子没有显含的系数。因此，固定后，只需耦合径向系数和有限个物理场分量，各模式可以独立求解。几何梯度引入的分量耦合已经由和处理，不妨碍这一分块。[5, §2]

以不可压缩向量扩散为例，未知量是正则性分量及压力。速度的二阶扩散落到基，压力梯度先落到基，再由提升一次；质量项则需要两个转换矩阵。散度约束单独保留在一次微分的输出基。于是，稀疏性不仅要求正确的微分系数，还要求对每一行方程跟踪其正则性指标和参数。[5, 式(28)–(30)]

对固定最高张量阶数，这种径向块的非零元数为，可交给现成稀疏线性代数程序处理。Coriolis 项会产生角向耦合，所以在这些时间演化实验中被放入显式部分；论文并未实现一个能对任意角向变系数保持相同分块结构的隐式求解器。球面 Rossby 波是另一个单独组装耦合的特征值例子。[5, §2、§4.1]

时间推进采用多步 IMEX。原文的更新形式为

本式是历史层的最大滞后指标，与空间分辨率无关：CNAB2取，SBDF4取，缺少的项以零系数处理。CNAB2为二阶，SBDF4为四阶；系数包含相应时间步长。论文中的比较使用常步长，以便分析时间误差。只要步长与线性算子不变，同一左端矩阵就可以反复使用其分解；这说明预计算和稀疏装配为何适合长时间积分。[5, §6.1, 式(53)–(54)]

完整算法中更昂贵的部分通常是变换。原文实际使用稠密矩阵乘法，径向与角向的工作量分别为

合起来为

另加经度 FFT。当与均为时，这是对个空间自由度执行的变换。因此，“微分只有一条带”和“整个三维求解过程准线性”是两个不同命题，后一个并非本文结论。[5, §2]

按变换矩阵的复用方式估算，固定张量阶数时，预存径向矩阵主要需要空间，角向矩阵需要空间；这是由实现结构推得的量级估算。高效 BLAS 和并行转置可以改善实际速度，但预计算、存储与数据通信仍是完整方法的一部分。

4. 特征值测试检验了哪些结构

球面 Rossby 波。 首先只在上求解

取，解析频率为

作者固定，每个场分量使用512个自旋球谐，验证自旋微分与乘法矩阵的配合。相应特征频率达到机器精度；这一结果与解能够由有限个相邻球谐组合表示有关，不能推广成任意谱特征值问题的全部离散特征值都能达到机器精度。[5, §4.1, 式(10)–(14)、Fig. 1]

球 Bessel 问题。 对

径向解析解为，由确定。作者使用512个径向多项式、，约前一半低径向波数特征值很准确，而接近截断端的特征值仍有较大误差。第100个特征函数同时显示正确的球心行为。这一测试把“基函数保证正则性”与“需要足够分辨率才能逼近振荡特征函数”清楚地区分开来。[5, §4.2, 式(15)–(17)、Figs. 1–2]

不可压缩向量扩散。 方程为

以磁向量势和标量势替代，同样的代数结构描述 Coulomb 规范下的磁扩散，其中。作者测试了五类边界：[5, §4.3, 式(18)–(27)]

边界类型	原文施加的主要物理约束
无滑移
不可穿透、无应力	，两个切向应力分量均为零
磁场势匹配	球内磁场与的调和外场匹配，边界关系依赖球谐次数
完全导体	法向磁场与切向电场为零；文中在所选规范下以切向及的零条件实施
伪真空	切向磁场为零，并与散度规范相容

势匹配是依赖全局角向展开的非局部条件；伪真空是局部切向零条件。两者在某些正则性分量中的公式相近，但不是同一个物理边界。

无应力边界尤其说明了任意阶张量工具的用途。先形成

再施加

例如不可穿透条件在正则性系数上是

因此不能简单地把“正则性分量”设成零并称为“径向速度为零”。应力已经经过一次微分，边界评价须使用基，并用二阶张量的提取所需坐标分量。[5, 式(23)–(24)、(31)、(33)–(34)]

5. 多项式如何改变有限分辨率结果

作者比较与。对二阶方程，前者把残差放在输出基的最高阶模式，等价于替换最后一条内部方程；后者在未知量中显式增加，并以的末列把残差转到输出基。涉及其他正则性分量时，相应使用。散度约束不添加这类残差，因而在截断系数意义下精确施加。[5, §4.3]

对、每个变量256个径向系数的扩散测试，矩阵维数分别为与，差别来自三个附加边界自由度。在所测情形中多得到约20个准确特征值，准确特征值数量约增加10%，完全导体边界的改善尤其明显。这是特定有限维残差选择的数值收益，不能提升为对任意 PDE、参数与分辨率都成立的最优性结论。[5, Fig. 3]

这里也不宜把误差变化归因于“相同离散系统改用了另一组基”。Part-I已经指出，只有在有限残差空间保持不变时，单纯换基才不改变精确解；本测试实质上改变了允许截断误差出现的位置和形状。

6. 边值测试：从磁场恢复向量势

为发电机计算构造初值，需要从给定的求

并满足外部势场匹配。三个旋度关系中存在依赖，单独求旋度方程还留下规范自由度。作者保留两个独立旋度关系，用散度零替换冗余方程，并将右端从提升到一次微分的输出基，组装可直接求解的稀疏边值系统。[5, §5, 式(41)–(46)]

测试磁场由径向多项式和少量球谐组成，解析向量势已知。在、时，径向、余纬、经度三个物理分量的相对网格最大误差分别为

这里的相对误差是各分量的最大绝对误差除以该分量解析解的最大绝对值。它联合检验了规范约束、边界与正反变换，而不只是检查一个稀疏微分矩阵。[5, 式(47)–(52)]

7. 三个非线性基准及其证据边界

第一个基准：边界驱动的旋转流。 方程为

球面施加给定速度，球内从静止开始，最终达到稳态。作者以

为诊断量，在、、时报告

这一约十位有效数字的结果与已有 Chebyshev 实现一致。由于终态是稳态，测试主要检验空间离散、非线性输运、Coriolis 项与速度边界，对时间离散细节相对不敏感。作者也承认，按相近自由度比较，另一套 Chebyshev 代码可能略准确；不能把该算例写成 Jacobi 基全面优于 Chebyshev 基的证据。[5, §6.3, 式(57)–(61)、Table 3]

这一算例还显示，球坐标速度分量在球心附近可以随方向变化，而笛卡尔速度仍光滑。判断物理场是否奇异，必须连同坐标基一起看，不能只检查单个球坐标分量的连续图像。[5, Figs. 4–5]

第二个基准：旋转对流。 原文采用

参数为、、、；速度不可穿透且无应力，温度在外边界固定。适当初值下形成的行波，空间结构持续旋转，动能则近似恒定。[5, §6.4, 式(62)–(65)]

以、SBDF4和计算，得到

虽然诊断量几乎不变，解本身仍随时间传播，所以时间误差不能省略。在相同空间分辨率下，二阶 CNAB2的动能误差按下降，四阶方法明显改善结果。因此，这个基准中对既有结果的精度改进有一部分来自时间推进；单凭动能的差异无法判定空间基底的优劣。[5, Table 4、Figs. 6–8]

第三个基准：旋转对流发电机。 这里增加磁向量势与 Lorentz 力：

参数为、、、、。速度不可穿透且无应力，温度固定，磁场与球外势场匹配；是 Coulomb 规范中的标量势。[5, §6.5, 式(66)–(70)]

系统最终呈现振荡发电机状态。磁能的定义包含特定归一化：

在欠分辨的计算中，仅改变便可能产生明显不同的长期能量演化，甚至出现磁场难以维持的情况。这说明低分辨率结果对截断残差敏感；一条曲线偶然接近某个参考值，不能代替空间与时间收敛检查。[5, Fig. 11、Appendix E]

作者同时检查了诊断方法。用“常数加单一正弦”拟合能量会受到高次谐波、采样及时间区间影响，因此转而记录长期极值，定义

磁能同理。这里上横线是极值中点，不是时间平均；是半极差，不是第一 Fourier 模态振幅。原文报告

这些数字来自提高径向和角向分辨率、减小 SBDF4步长，以及对长期振荡极值的检查。各指标受到的截断和瞬态限制并不完全相同，原文还指出部分量若要进一步提高精度，需要继续消除长时间瞬态的影响；不宜将全部指标统称为具有统一的七位精度。[5, 式(71)–(81)、Tables 5–6]

8. 并行结果与本篇真正增加的贡献

在的对流强扩展测试中，原文给出：

核数	本文球内实现每步时间	作为参照的 Rayleigh 每步时间
64	6.22 s	2.42 s
1024	0.53 s	0.19 s

从64核增至1024核，本文相对加速约11.7倍，理想值为16倍；实际每步时间约为参照代码的三倍。不过，Rayleigh 在这里计算球壳，可采用径向 DCT；本文计算包含球心的完整球，径向也使用矩阵变换。原文还指出连续分配区块会造成负载不均衡。因此，这个表格说明实现具有良好强扩展，同时暴露了变换与负载分配的成本；它不是同域、同算子、同自由度的严格速度竞赛。[5, §6.2、Table 2]

Part-II最重要的贡献，是把 Part-I 的几何结构落实为一条可以逐项检验的计算流程：预计算和变换处理动态范围；稀疏系数装配跟踪正则性与基参数；边界检验覆盖向量与应力；非线性基准进一步分离空间、时间、截断残差和诊断误差。作者提供的证据支持这套方法能够可靠处理复杂球内张量 PDE，但不支持“所有特征值都准确”“所有边界选择都同样好”或“比所有现有代码更快”这样的概括。

两条路线的深入比较

这两条路线共享许多基础：都利用圆盘或球的旋转几何，都需要角向 Fourier 或球谐展开，都借助正交多项式递推与带状矩阵，也都把几何上复杂的操作改写成一维结构化运算。它们的主要区别，在于首先优化的对象不同。

DFS–BMC路线首先构造一个可自适应计算的函数表示：

此处主要的压缩对象是不同坐标之间的相关性。一个给定函数若只需要很小的，求值、积分和许多函数运算都能复用少量一维对象。

Jacobi路线首先构造几何相容的有限维空间及空间之间的算子：

此处主要的结构是微分所要求的局部正则性与参数转换。它不假定函数在径向与角向之间具有很低的分离秩，而是尽量使任何合法谱系数都能被统一微分。

因此，低分离秩与物理张量阶数不能混用。前者描述坐标分离表示有多少项，后者描述物理量有多少个张量指标。一个向量场可以具有很高的分离秩；一个标量场同样可以具有很高的分离秩。文献[3]–[5]所说的 tensor calculus 主要是物理张量微积分，并不表示已经对三维系数数组作了低秩张量分解。

比较层面	DFS–BMC与低秩路线：[1]–[2]	正则性适配的 Jacobi 路线：[3]–[5]
首要目标	自适应表示函数，并支持广泛函数运算	统一标量、向量、张量的稀疏谱微积分与 PDE 装配
几何处理	双重覆盖、滑移反射、BMC对称性及极点或圆心值约束	端点幂、自旋、正则性分量和相应 Jacobi 参数
主要压缩或稀疏来源	坐标分离秩小、FFT/DCT、特定算子的带状结构	微分升降只改变少量指标，基转换与简单乘法带状
未知函数的组织	少量自适应行列切片及其一维谱表示	按角阶、正则性、径向次数组织的完整系数族
向量与张量	论文实现用笛卡尔分量支持向量运算	从旋转或自旋基出发，给出任意阶张量的统一规则
PDE证据的主体	球面和圆盘上的快速标量 Poisson 求解	圆盘线性问题，以及球内特征值、边值和非线性演化
主要成本风险	分离秩增长、主元搜索及重压缩、操作后正则性控制	稠密变换和预计算、存储、通信及较复杂的分量组织

表中后两列概括的是这五篇论文的实际方法与测试范围，不能据此断言某条路线在后续研究中无法支持其他问题。

下面对这两条路线进行部分更细致的比较：

• BMC 对称性与完整正则性并不等价：这是比较中最需要仔细区分的一点。BMC保证参数域上的两份数据对应同一个物理点，BMC-I或BMC-II进一步保证极点或圆心函数值一致。但是，函数值一致与全部笛卡尔导数具有正确极限，是不同层次的要求。在圆盘上，对解析标量，角向 Fourier 系数必须具有
  
  的结构。它同时包含奇偶性、低阶 Taylor 系数消失以及所有更高阶相容性。相反，BMC-II本身主要给出径向奇偶关系与原点常值关系。下面的反例是本综述为说明这一差别给出的推导，不是论文中的数值算例。取
  
  它在加倍矩形上解析，并且
  
  所以它满足BMC-II，奇偶性正确，非零角频也在原点取零。可是回到笛卡尔坐标后，
  
  这个函数在原点连续，却不可微。因为、，若可微，则梯度只能是；但，与沿对角线应有的一阶近似矛盾。其分量只有阶消失，而完整正则性要求。由此得到的准确结论是：保BMC的低秩截断，并不单凭对称性就保证每一个有限截断都属于笛卡尔光滑函数空间。文献[2, §2.4、§4.4、§5]中从BMC-II结构及径向可除性推断物理圆盘上连续可微的论证，还需要角模态与径向消失阶之间的进一步相容条件。仅有奇偶性、原点常值以及一次除以的可除性，并不足以保证笛卡尔意义下的连续可微。这并不否定DFS对真正光滑输入函数的高精度逼近。若输入本来就是光滑的笛卡尔函数，它的精确拉回当然具有正确的相容关系；[1]–[2]还专门设计了对称性维护、极点或圆心的低阶处理，以及系数空间中的除或除操作。需要区分的是：对一个光滑目标获得小函数误差，与证明任意截断表示在物理坐标下具有所有阶光滑性，不能互相代替。Jacobi路线在这一层具有更直接的空间保证。圆盘基自带；球面自旋球谐自带正确的半角端点幂；球内经变换后的各分量自带。只要保留合法模式并按既定规则转换分量，每个有限展开都保持相应正则性。这尤其有利于反复求导和张量分量耦合[3, §2–§3；4, §2–§4]。
• 函数逼近 — 低秩自适应与几何适配：对于许多待计算函数，径向与角向变化并不强烈耦合。例如轴对称标量在角向只需一个模式，分离秩很小。DFS–BMC路线能够从可调用函数中自适应发现这一点，并分别决定一维分辨率与秩。它的优势不是必须先把整个二维张量积网格做满，而是能在适当条件下只处理少量关键切片[1, §3；2, §3]。这一收益取决于函数和所用坐标。光滑性并不保证低分离秩；沿不同方向传播的振荡、多个尺度的角径耦合以及后续非线性乘积都可能增加。加法和乘法也会先扩大表示，再需要重压缩。连续 GE 的误差收敛理论有其解析性和主元选择假设，实际自适应采样并不是对任意黑盒函数的无条件误差证书。Jacobi路线节省的是另一类工作：避免用大量不满足几何相容性的模式去表示应当高阶消失的角分量。角阶越高，合法径向空间越受到或限制，按总次数作三角截断可以相应减少径向自由度。它不需要依赖数值相消来重建这些端点幂，但也不能像低秩表示那样自动利用任意坐标之间的可分离性[3, §2、§5.3；5, §2]。两者都依赖解的实际正则性才能获得谱精度。Jacobi基满足解析场的几何结构，不表示带真实物理奇性、间断或未分辨边界层的解也必然指数收敛；DFS延拓同样不能把一个本来不光滑的物理解变成光滑解。这里必须把坐标造成的表面奇点与 PDE 解自身的真实奇性区分开。
• 微分与稳定性：在笛卡尔光滑场上，球坐标表达中的大项往往需要相消。例如
  
  若在网格上分别计算导数与项，靠近原点时可能失去有效数字。Jacobi路线把这种相消事先完成在基函数恒等式中，得到单带；导数输出进入正确的新正则性和新参数空间。这是它在高角阶与重复张量微分时的实质优势。DFS–BMC路线也没有简单地在极点网格上直接除零。论文对球面切向导数和圆盘笛卡尔导数使用系数空间的带状乘法矩阵求解除或除，并配合奇偶性及零值条件。这比直接对退化坐标公式逐点计算更可靠；不过，这些低阶除法结构与所有高阶模式相容性仍是不同问题[1, pp. C416–C420；2, §4]。另一方面，稀疏不等于条件数小。Jacobi微分系数随次数增长，边界评价行并不局限于少数系数，完整 PDE 也可能具有病态参数、近零特征值或非正规性。128位预计算、残差选择以及高频特征值检查，正说明仅有漂亮的稀疏矩阵还不足以结束稳定性分析[4, §5；5, §3–§4]。因此，较客观的表述是：Jacobi路线在离散空间层面更直接地避免坐标正则性失配；DFS–BMC路线以结构化延拓和专门系数操作获得有效算法。某一个 PDE 在给定精度下哪种实现更稳定，仍应检查相应离散算子、边界与误差，而不是只依据“Fourier”或“Jacobi”这个基底名称判断。
• 求解成本：先看二维函数表示。若两个方向分别需要个系数，DFS–BMC的秩表示需要存储，其一维因子的FFT或DCT工作量随这些因子的长度增长。当时，这会明显少于完整二维数组。若接近其中一个方向的分辨率，优势便相应减弱，主元搜索和重压缩也成为额外工作[1, §3；2, §3]。Jacobi圆盘算法对每个角模态执行稠密径向变换；径向和角向分辨率同阶时，这部分通常为。但对已转到系数空间、系数结构合适的线性 PDE，每个角模态的径向算子很稀疏，求解可以接近自由度数的线性成本。前者是变换费用，后者是求解费用，不能用后者覆盖整个流程[3, §5.3、Example 4]。对三维球内实现，文献[5]的变换作用在个自由度上。这是其真实成本边界；但不能拿它直接与文献[1]二维球面上的系数求解比较，从而宣布球内方法更慢。三维问题本身多了一个径向变量。在真正共同的球面标量 Poisson 问题上，两条路线的取舍更清楚。由文献[4]的球谐关系直接推出
  
  兼容条件是；指定零均值时取。所以球谐系数中的 Poisson 求解本身是对角的，主要工作在正反变换。Spherefun则用双 Fourier 表示和带状纬向方程完成同一任务，兼容其自适应函数对象，并可利用 FFT。这个比较说明，“最稀疏的系数求解”与“最快或最方便的完整函数工作流”未必由同一个表示同时实现。圆盘上也有类似区别。Diskfun先把齐次边界因子提出，再利用 Sylvester 方程及低秩ADI；若右端与解具有较小分离秩，能够节省存储与工作。Jacobi方法按角模态使用几何适配的径向带状系统，不要求右端低分离秩。ADI要计入迭代数、非正规谱信息与重压缩，Jacobi方法要计入变换与预计算；两者都不能只报最有利的内部求解步骤[2, §5；3, §5]。
• 向量、约束和边界条件：前两篇已经支持向量计算，不能将其概括为只能处理标量。其做法主要是分别表示笛卡尔分量，再构造梯度、散度、旋度和相关操作。这种方式与用户熟悉的笛卡尔物理量相近，且避免将球坐标基向量的表面奇异性暴露给每个分量[1, §4.1；2, §4]。后三篇的突出优势，是把分量变换与协变微分系统化。圆盘上的、球面上的自旋、球内的正则性都被明确跟踪。速度梯度成为二阶张量，应力再由对称化得到，散度或旋度通过既定收缩关系组合。这种组织更便于保留压力、速度或磁向量势等原始变量，并逐项添加新的张量方程[3, §3；4, §3；5, §4.3]。但是，通用张量表示也有代价。若不利用对称性，三维阶张量有个分量；“任意阶”不表示成本与张量阶数无关。无散向量的 toroidal–poloidal 标量势表示能够减少独立变量，专门设计的标量势代码在某些问题上仍可能更经济。文献[5]的原始变量方法更便于统一表达方程与边界，而不是保证对每个特殊系统都使用最少未知量。边界方面，Diskfun的因子以及 Jacobi 的基重组，都利用了“乘一个边界消失因子”这一思想；两条路线在这里并非完全分离。Jacobi文章进一步详细比较了残差的选取，并通过向量、应力和磁边界给出系统验证。其困难在于边界迹必须使用正确的分量基与导数后的参数，不能把某个正则性分量直接当成法向或切向分量[2, §5；3, §5.2；4, §5；5, §4.3]。
• 非线性 PDE 与时间演化：就这五篇论文的实际完成程度而言，Jacobi路线对非线性张量 PDE 的证据更完整。Part-II给出了网格乘积、MPI变换、IMEX、不可压缩约束及磁规范的组合，并以三个完整基准检验长期演化[5, §6]。前两篇的主体则是函数计算和快速 Poisson 求解。它们提供了后续构建演化算法所需的许多工具，但不能把函数对象支持乘法与微分直接等同于“文中已经验证了任意非线性球面 PDE 的自适应低秩积分器”。特别是非线性会耦合原本分开的坐标结构，分离秩是否保持小，需要对具体问题另行分析。Jacobi路线也没有绕开全部困难。网格乘积需要处理混叠，显式部分仍受到时间尺度约束，原文三维实现不能将任意角向耦合都放进原有的隐式分块。对快速旋转系统，把Coriolis项显式处理可能影响允许步长，这是从其显隐式划分得出的算法层面考虑，而非论文给出的普遍性能定理。Part-II尤其有价值的一点，是说明精度差异未必来自空间基底：稳态流体例子主要检验空间与边界；传播对流例子突出时间误差；发电机例子同时暴露欠分辨残差和诊断方式的影响。客观评价一种方法，应当先控制这些因素，再比较基底本身。
• 场景适配，以及可能的互补方式: 若主要任务是从可调用的光滑函数出发，自动决定分辨率，反复求值、积分、可视化并执行函数代数，而函数又具有较低坐标分离秩，那么前两篇的设计更直接。用户得到的是能够在许多操作中复用的函数对象，不需要先为某一PDE手工固定全部模式数。若主要任务是长期反复求解包含向量、压力、应力、磁场或更高阶张量的圆盘、球面或球内方程，尤其需要严格组织坐标奇点正则性和微分后的分量空间，那么后三篇的框架更系统。它的额外实现与预计算成本可以在多次算子作用和时间推进中摊销，但仍应预计稠密变换与通信的负担。对于规则几何上的常系数标量椭圆问题，两条路线都有十分高效的专门结构，不能仅凭上述一般偏好决定。应具体比较输入从何而来、是否已在某种谱系数中、是否需要自适应函数表示、右端与解的分离秩、边界类型，以及重复求解次数。两条路线还存在一个合理的互补方向：利用自适应函数采样发现尺度与分辨率，再投影到正则性适配的 Jacobi或自旋球谐空间中执行 PDE 运算；也可以在几何合法的系数表示内部研究额外的低秩压缩。不过，这需要控制投影误差、保留各角模态的端点幂，并核算两套表示之间的转换费用。

综合这些证据，两条路线各自最鲜明的贡献可以表述为：前两篇把圆盘和球面的函数计算组织成可自适应压缩、可复用的一维计算；后三篇把退化坐标中的物理张量微积分组织成保持正则性的稀疏谱运算。前者的实际收益主要取决于分离秩和函数工作流，后者的实际收益主要取决于几何相容性、算子重复使用及张量复杂程度。公平的比较应以同一物理问题、相同误差标准和完整计算流程为基础。