Volterra Integral Equations | ALL YOU NEED IS LOVE

考虑如下形式的Volterra积分方程：

其中，是三角形区域上的一个给定核函数。下面介绍几种常用来求解此种问题的数值方法。

Collocation Method

REFERENCES
Shen, Jie, Tao Tang, and Li-Lian Wang. Spectral methods: algorithms, analysis and applications. (2011) Section 5.1

配点法的思路非常直接：选定一组配点，要求找到一个多项式使得在这些配点处满足方程，即

然而由于方程中的积分是变上限的，这样的形式不利于构建离散的方程，为了对于近似点值建立线性系统进行求解，我们需要先使用变量代换来将上式的变上限积分转化为固定的上的积分。令

其中而，于是原方程变为

对应的配点方程为

下一步我们需要建立插值矩阵，使用配点网格上的值进行插值来近似，从而将上式中的积分转化为一个矩阵-向量乘积。设

是以配点网格上的Lagrange基函数为基底的插值多项式，则

接下来将上面的表达式代入配点方程中得到

我们可以把这一表达式写为矩阵形式：

其中是第个标准基向量，，。将对应的方程合起来，我们就得到了一个线性系统

其中

一般需要通过数值积分来计算，常用的数值积分方法包括高斯-勒让德积分和Clenshaw-Curtis积分等，因此矩阵的构建复杂度依赖于核函数的光滑性以及数值积分方法的效率，通常复杂度为，其中是计算积分使用的点的数量(相比于一般无法忽略)。

通常实际计算中为了兼顾计算效率与数值稳定性会将配点选为Gauss点，此时基函数可被相应的正交多项式表示，例如当使用Legendre-Gauss配点时，可以将上述积分中的Lagrange基函数表示成Legendre多项式展开的形式：

其中当使用Legendre-Gauss与Legendre-Gauss-Radau配点时，，而当使用Legendre-Gauss-Lobatto配点时，，对于。

最后我们可以通过求解上述线性系统来得到近似解。由于矩阵通常是满矩阵，这类方法的总体复杂度为，其中是配点的数量。

Chebfun目前使用此类方法来求解Volterra积分方程，具体实现细节可以参考Chebfun的文档和源代码。

Gutleb-Olver’s Method

REFERENCES
Gutleb, Timon S., and Sheehan Olver. A sparse spectral method for Volterra integral equations using orthogonal polynomials on the triangle. SINUM (2020).

Gutleb-Olver的方法的核心思想是在将Volterra积分方程中的核函数在三角形区域上展开成一组适当的正交多项式的级数的前提上，首先将需要积分的原本定义在一维区域上的函数表示成一个同样定义在三角形区域上的二维函数，要求

再使用同样的三角区域上的正交基将展开，接着计算三角区域上的函数乘积并关于积分，由此将原方程中的Volterra积分转化为一个矩阵-向量乘积，最后组装线性系统并进行求解。为了完成这一套操作，Gutleb-Olver引入了以下几个算子：

延拓算子：将定义在上的函数延拓为定义在三角形区域上的函数，要求。
乘积算子：将定义在三角形区域上的函数乘以核函数得到一个新的函数。
积分算子：对定义在三角形区域上的函数关于第二个变量进行积分得到一个新的函数。这里的是一个权重矩阵，后面我们将看到它实际上是对应的乘积算子，是关于变量的积分算子。

通过上述三个算子，我们可以将原方程中的积分项转化为一个矩阵-向量乘积，从而得到一个线性系统来求解近似解。

为了表示定义在三角形区域上的函数，Gutleb-Olver使用了Koornwinder多项式作为正交基：

其中是Jacobi多项式。对于任意定义在三角形区域上的函数，我们可以将其展开成Koornwinder多项式的级数：

通过定义

我们可以将上述函数展开写为
$，$
此处包含了所有全阶数(关于的阶数与关于的阶数之和)为的Koornwinder多项式，是对应的系数向量，而是包含了所有阶数的Koornwinder多项式系数的列向量。我们使用三角形角标表示该向量是三角形区域上展开的系数向量。

对于一维区间上的函数，Gutelb-Olver使用指标为的Jacobi多项式作为正交基，即最终的解函数的展开形式为

其中是将重新放缩到下的对应正交基函数。这一选择有多重原因，我们后面将看到这一基函数所带来的种种便利。Gutleb-Olver在这篇文章中的主要工作是分析出这组基下的算子、和的无穷维矩阵表示，发现在基下恰好是一个对角矩阵，并基于这一观察给出一种基于算子版本的Clenshaw算法快速构造Volterra积分算子的矩阵表示

的方法，于是以为核函数的Volterra积分算子在这组基下有以下无限维矩阵表示：

Volterra Integral Operator

为了构造Volterra积分算子，我们需要一些准备工作。首先定义二元函数在变量下的积分算子：设是定义在三角形区域上的函数，则对于每个固定的，我们可以将看成是变量上的一个函数，因此我们可以定义如下的积分算子：

另一方面，为了让一维的输入函数可以与二维函数进行乘积，我们需要将延拓为一个定义在三角形区域上的函数，要求，因此我们可以定义如下的延拓算子：

把这两个算子结合起来我们可以用，和表示算子对应的系数矩阵：

在上面的分析基础上，Gutleb-Olver方法最重要的观察是注意到

上述恒等式表明在特别挑选的基下，有以下事实：

其中是在下的三对角乘积矩阵，是一个对角矩阵。

现在令，对比公式和，我们可以发现实际上是已知的对角矩阵，这在后续的Volterra矩阵构造中非常重要：假设核函数可以在单项式基下被展开为

则对应的Volterra积分算子可以使用上面定义的辅助算子们来表示为

而这其中的分别是三角基下的和的分块三对角乘积算子：

具体表示为

它们满足以下交换性质：

此处的是基下的的三对角乘积算子：

利用上述交换性质，我们可以将Volterra积分算子表示为

不过由于核函数是由三角基表示的，对应的Volterra算子更适合写成

的形式，无法直接将三角基下的乘积算子中的和替换成，因此Gutleb-Olver基于三角基函数所满足的三项递推关系，使用算子版本的Clenshaw算法来进行快速计算。

Clenshaw Algorithm

Clenshaw算法是一种用于高效计算正交多项式级数的方法，特别适用于当正交多项式满足三项递推关系时的情况，经典一维Clenwshaw算法的讨论参见Clenshaw算法。对于二维的三角基函数，由于它内部是两个耦合的Jacobi多项式的乘积，所以同样满足类似的三项递推关系，因此和一维类似，对于任意三角区域内的点对都存在一个分块上Hessenberg矩阵满足

于是如果核函数，则

现在把上述点值等式提升为算子版本，我们可以得到

上述表示已经足够紧凑，但是还未利用上是对角矩阵的核心优势，下面我们记

并令

由于作为三角基函数具有三项递推关系，所以对应的也满足对应的关系，因此我们同样可以构造出一个分块上Hessenberg矩阵满足

这里右端项中是而非是因为当，即时，对应的。此处形如

其中表示一维系数空间上的恒等算子，是阶单位阵，与定义为

现在我们有

进而当时，下式成立：

通过代数变换，我们可以将这一等式改成更符合Clenshaw算法的等价形式，即

由于是一个分块上Hessenberg矩阵，基于上面公式可以高效得到的矩阵表示，进而得到Volterra算子。

注：在普通的Clenshaw算法中，是数值矩阵：

但在算子版本的Clenshaw中，第层未知量是由矩阵块组成的向量：

其中每个都是一维 Jacobi 系数空间上的矩阵算子，因此，表示用的数值系数混合这些算子块：

这里的表示不改变每个矩阵块本身。符号表示左乘一维 Jacobi 乘法矩阵：

因此

也就是说，先用混合 operator blocks，再对每个结果左乘。

类似地，表示右乘：

因此

所以实际上我们隐式地在的夹心结构中使用了交换性质：

求解方程

Gutleb-Olver 方法首先构造的是积分区间为到的核心算子，它对应的积分表达式为

完整的到的积分矩阵还需要乘上权重，即

若要求积分区间为到，可以使用反射变换：令新的变量为，则

因此，可以先在反射变量中构造积分限为到的Volterra 矩阵，将核函数改为

于是反射坐标下的矩阵为

如果要把结果重新表示回原来的坐标，还需要在输出端作用反射算子：对于上的Jacobi多项式，反射对应交换两个 Jacobi 指标：

因此，反射算子把函数变成，同时把系数从空间转换到空间，并带上符号因子，因此通过允许基变换，实际上是两组基之间的一个对角系数变换矩阵。Gutleb-Olver 用这个算子把的积分形式转换为的 Volterra 积分形式：

等价地，在算子版Clenshaw 递推中，可以不显式改写核函数，而是把所有对应变量的左乘算子替换为

而变量对应的右乘算子保持不变：

也就是说，算子对应的矩阵的构造可以理解为：在 Clenshaw 递推中把核的第一个变量反射为，然后再乘上权重，最后根据需要作用输出反射算子。

对于第二类方程，为了让积分项和恒等项仍然处在同一个系数空间中，通常还需要使用降阶和升阶算子进行基底匹配。升阶算子用来提高 Jacobi 指标。例如

以及

它们只改变系数表示，不改变函数本身：

降阶算子一般不是无权降阶，而是带权降阶，也就是说，降低指标时通常要同时乘上或，否则对应的矩阵不再稀疏。论文中给出的典型形式是

以及

所以的作用不是单纯把转成低参数 Jacobi 基，而是把乘权重后的函数表示在低参数空间中。

在第二类 Volterra 方程中，反射后若直接相加，会导致恒等项和积分项处在不同 Jacobi 空间中。因此 Gutleb-Olver 使用组合

从而把积分项重新放回系数空间中。这里吸收了权重并降到空间，做反射，最后升回空间。

综上，第二类方程

在Gutleb-Olver方法下的对应离散系统为

由于的矩阵表示可以通过Clenshaw算法高效计算，而最终的系数矩阵是严格带状的，构造复杂度为，求解复杂度为，其中是截断维数，为核函数的深度，对应最终系统的带宽。

选取基底的理由

Gutleb-Olver 方法并不是数学上只能使用，而是在其三角 Jacobi 基框架下，是保持稀疏结构的自然选择。他们在三角区域

上使用标准三角 Jacobi 基。其基函数可写为

对积分时，有

由于是上的Legendre 多项式，

所以只有的模态保留下来，从而

因此，的自然输出形式是，这说明权重因子为，而去掉权重后的只需要在每个总阶数固定的块中选取第一个分量。更重要的是，若从一维空间延拓到三角基，则

是对角矩阵，这一对角性是 Gutleb-Olver 快速构造的核心。它使得算子版本的Clenshaw递推可以直接从出发，而不需要显式构造完整的三角乘法算子。此外，对于第二类 Volterra 方程，未知量同时出现在恒等项和积分项中，如果输入基和输出基不同，即公式()中的非零，那么恒等项就需要额外的连接矩阵。一般情况下，这种转换矩阵不是固定带宽稀疏矩阵而是稠密的上三角矩阵，从而将破坏整体带状结构。

总地来说，选择有以下几点原因：

在三角基下，沿积分后自然得到；
具有简单的块选择结构；
是对角矩阵，算子层面的Clenshaw算法可以利用这个对角结构快速构造；
第二类方程中恒等项和积分项可以保持在同一个系数空间，由此得到的离散系统保持稀疏/带状结构。

换一种角度来想，如果希望使用其他指标为基底，立刻会发现为了保持系统的稀疏性，有以下限制存在：

为了让()成立以保证是对角矩阵，必须令；
为了避免使用稠密的基转换矩阵，必须令；

这样就唯一确定了输入基只能为。