第二类Volterra积分方程的广义多步配置法

Generalized Multistep Collocation Methods for the Second-Kind Volterra Integral Equations

摘要: Volterra积分方程作为一种重要的数学模型，被广泛应用于多个科学研究领域。针对第二类Volterra积分方程的数值解，本文在经典多步配置法的基础上，结合边值方法的思想，研究出一种广义多步配置方法。该方法利用Lagrange插值公式，以不同的节点作为插值节点，将原方程离散成为一个线性方程组。通过实验，本文验证了该方法在求解第一类Volterra积分方程的有效性，并且可以达到较高的收敛阶。Volterra积分方程在实际科学领域中有广泛的应用，因此对其数值求解方法的研究具有重要的理论意义和应用价值。本文所提出的广义多步配置方法为求解第二类Volterra积分方程提供了一种高效且可靠的数值逼近方案，为相关领域的研究和应用提供了新的思路和方法。

Abstract: As an important mathematical model, Volterra integral equation is widely used in many scientific research fields. Aiming at the numerical solution of the second kind of Volterra integral equation, this paper studies a generalized multi-step collocation method based on the classical multi-step collocation method and the idea of boundary value method. This method uses the Lagrange inter-polation formula to discretize the original equation into a linear equation set with different nodes as interpolation nodes. Through experiments, this paper verifies the effectiveness of the method in solving the first kind of Volterra integral equation, and can achieve higher convergence order. The Volterra integral equation has a wide range of applications in the field of practical science, so the research on its numerical solution method has important theoretical significance and application value. The generalized multi-step collocation method proposed in this paper provides an efficient and reliable numerical approximation scheme for solving the second kind of Volterra integral equations, and provides new ideas and methods for research and application in related fields.

1. 引言

在许多实际问题中，如物理学、生物学、化学、经济学等多个领域，Volterra积分方程都有着十分广泛的应用。常见的热传导模型、流行病扩散、声学散射问题等都可以用Volterra积分方程来模拟 [1] [2] [3] [4] 。然而，对于一般的Volterra积分方程，其真实解往往难以通过解析方法得到，因此，如何通过数值方法求解Volterra积分方程的近似解成为了一个重要的研究课题。近年来，许多学者对数值求解Volterra积分方程进行了深入研究，提出了一系列有效的数值方法，例如线性多步法、Runge-Kutta方法、谱方法、配置法等。这些方法的应用范围广泛，可以有效地逼近Volterra积分方程的解，为解决实际问题提供了有力的工具。研究Volterra积分方程的数值求解方法具有重要的理论和应用价值。它可以为解决实际问题提供有效的数学模型和算法，促进各领域的理论研究和实际应用的发展。

配置法是一种有效的数值求解Volterra积分方程的方法，它通过选择一组适当的配置点，利用这些配置点的特性来逼近Volterra积分方程的解。配置法不仅具有较好的数值精度，而且计算量相对较小，因此在求解Volterra积分方程时具有很大的优势。多步配置法是将配置法与其他数值方法相结合的一种高效数值求解方法。由于多步配置法具有收敛性好、易于构造等优点，因此广泛受到国内外学者的重视。在多步配置法中，通常选择适当的配置点，并根据这些配置点的特性来构造一个多项式，从而得到Volterra积分方程的解的近似值。通过选择不同的配置点和多项式，可以获得不同的多步配置法。

在求解Volterra积分方程时，多步配置法的计算精度和计算量都可以达到很高的水平。与其他数值方法相比，多步配置法的收敛速度通常更快，因此在求解大型Volterra积分方程时具有很大的优势。此外，多步配置法的适用范围也较为广泛，可以用于求解各种类型的Volterra积分方程，包括线性、非线性、离散、连续等类型的Volterra积分方程。因此，多步配置法已成为求解Volterra积分方程的一种重要数值方法。

1983年，Wolkenfelt运用直接求积法求解第二类Volterra积分方程，并讨论其数值稳定性 [5] 。1984年，Brunner利用配置法和迭代配置法来离散并求解第二类Volterra积分方程，并给出了这两种方法的收敛性分析 [6] 。1994年，韩国强利用迭代配置法求解了第二类Volterra积分方程，通过在节点处利用配置解的渐进展开进行Rihcardson外推，从而提高了配置解的精度 [7] 。2012年，Liang和Brunner研究了第一类Volterra积分方程的配置方法，并对其配置解的局部超收敛性进行了分析 [8] 。2013年，Xiang和Brunner提出了高效的Filon配置方法用于求解Volterra积分方程 [9] 。2016年，Zhang，Liang和Brunner提出了第二类自卷积Volterra积分方程的配置方法，证明了数值解的存在唯一性，并证明了该方法是收敛的 [10] 。

经典的配置方法在求解Volterra积分方程时具有较好的数值精度和计算效率。然而，当需要达到高精度时，配置点的个数会相应增加，从而导致计算量的增大。为了克服这一困难，许多学者开始考虑在不增加配置点的情况下提高算法收敛阶的数值方法。

2009年，Conte和Paternoster利用单步方法得到的近似值构造了多步配置方法。该方法能够在不增加配置点的情况下提高了多步配置法的收敛速度和精度，同时对该方法的收敛性及其线性稳定性进行了分析 [11] 。2012年，Fazeli等讨论了多步配置方法求解第二类Volterra积分方程的稳定性，给出了方程的超隐式配置解，且这类数值格式具有较广的稳定域 [12] 。2015年，Fazeli和Hojjati介绍了Volterra积分微分方程的超隐式多步配置方法 [13] 。2017年，Ma和Xiang基于特殊的多步配置方法，利用未计算的近似值，提出求解Volterra积分方程的配置边值方法，并分析了其线性稳定性 [14] 。2018年，Zhang和Liang将多步配置法应用于第一类Volterra积分方程，证明了多步配置解的存在唯一性，并分析了该方法的收敛条件，给出了相应的收敛阶 [15] 。2019年，Zhao和Long等人利用光滑变换，将弱奇异Volterra积分方程转化为具有更好正则性的弱奇异Volterra积分方程，变换后的积分方程的解充分光滑。然后利用已计算的近似值和当前以及下一个子区间的配置点构造了超隐式多步配置法，有效地减少了计算量和内存需求。此外，该方法还具有良好的数值稳定性和收敛性 [16] 。2020年，Chen，Liu和Ma [17] 给出了第二类高振荡Volterra积分方程的一般多步方法，并分析了该方法的收敛性和稳定性。2022年，Zhao和Fan等针对带有高振荡核的线性Volterra积分方程提出配置方法，采用了Filon型方法对配置方程中的振荡积分进行离散，并研究了该方法的收敛性 [18] 。

本文针对第二类Volterra积分方程，提出了一种广义多步配置方法，旨在解决大规模的复杂问题。该方法结合了配置法和多步法的优点，通过选择适当的配置点构造广义多部配置法逼近第二类Volterra积分方程的解。

2. 算法构造

在这一节中，我们将构造第二类Volterra积分方程的广义多步配置法( $SGM{C}_{k_1-m-k_2}$ )。考虑如下第二类Volterra积分方程

$u\left(t\right)=f\left(t\right)+{\displaystyle {\int }_0^{t}K}\left(t,s\right)u\left(s\right)\text{d}s,t\in \left[0,T\right].$ (1)

其中， $u\left(s\right)$ 为未知函数， $f\left(t\right)$ 和 $K\left(t,s\right)$ 为已知函数，且 $f\left(t\right)\in C\left(\left[0,T\right]\right)$ ， $K\left(t,s\right)\in C\left(\Delta \right)\left(\Delta :=\left\{\left(t,s\right):0\le s\le t\le T\right\}\right)$ ，并有 $f\left(0\right)=0$ ， $\left|K\left(t,t\right)\right|>0$ 。

对方程(1)的区间 $\left[0,T\right]$ 进行均匀网格划分，利用节点 $t_n=nh,n=0,1,\cdots N$ ，将区间 $\left[0,T\right]$ 分割为N个子区间。选择均匀网格中每个小区间 $\left[t_n,t_{n+1}\right]$ 上的端点，并在区间 $\left[t_n,t_{n+1}\right]$ 前取 $k_1$ 个节点，区间 $\left[t_n,t_{n+1}\right]$ 内取m个点，区间 $\left[t_n,t_{n+1}\right]$ 后取 $k_2$ 个节点，利用这些节点作为插值节点构造Lagrange插值函数去近似未知函数。具体操作步骤如下：

首先对区间 $\left[0,T\right]$ 进行均匀网格划分

$I_h=\left\{t_n:t_n=nh,n=0,\cdots ,N,h\ge 0,Nh=T\right\}.$

这样就得到N个子区间，其中h为步长。定义每个子区间 $\left[t_n,t_{n+1}\right]$ 内的配置点为

$t_{n,j}=t_n+c_{j}h,j=1,\cdots ,m\left(0<c_1<\cdots <c_m<1\right).$

方程(1)可以改写为

$f\left(t\right)=F_n\left(t\right)+{\Phi }_n\left(t\right),t\in \left[0,T\right],$ (2)

其中

$F_n\left(t\right)={\displaystyle {\int }_0^{t_n}K\left(t,s\right)u\left(s\right)\text{d}s}$

被称为滞后函数。

${\Phi }_n\left(t\right)={\displaystyle {\int }_{t_n}^{t}K\left(t,s\right)u\left(s\right)\text{d}s}$

被称为增量函数。

定义在节点 $\left\{t_k,c_j|k=-k_1,\cdots k_2+1,j=1,\cdots ,m\right\}$ 处的Lagrange插值基函数分别为 ${\varphi }_k^{k_1,k_2}\left(s\right)$ 和 ${\psi }_j(\; s\; )$

${\varphi }_k^{k_1,k_2}\left(s\right)={\displaystyle \underset{\begin{array}{l}i=-k_1,\\ i\ne k\end{array}}{\overset{k_2+1}{\prod }}\frac{s-i}{k-i}}{\displaystyle \underset{i=1}{\overset{m}{\prod }}\frac{s-c_i}{k-c_i}},k=-k_1,\cdots ,k_2+1,$

${\psi }_j\left(s\right)={\displaystyle \underset{i=-k_1}{\overset{k_2+1}{\prod }}\frac{s-i}{c_j-i}}{\displaystyle \underset{\begin{array}{l}i=1,\\ i\ne j\end{array}}{\overset{m}{\prod }}\frac{s-c_i}{c_j-c_i}},j=1,\cdots ,m.$

接下来构造在区间 $\left[t_n,t_{n+1}\right]$ 上的配置多项式，有如下形式：

当 $0\le n\le k_1-1$ 时，在区间 $\left[t_n,t_{n+1}\right]$ 上的配置多项式为：

$u_h\left(t_n+sh\right)={\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1,k_2}\left(s\right)y_{k_1+k}}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\psi }_j\left(s\right)U_{n,j}},s\in \left[0,1\right].$

由于当 $0\le n\le k_1-1$ 时向前取到的配置点个数不足，无法取到全部的配置点，为了确保能够取得所需数量的配置点，我们向后借用后面区间端点为配置点。这样得到的配置点分别 $y_0,y_1,\cdots ,y_{k_1+k_2+1}$ 和 $U_{n,1},U_{n,2},\cdots ,U_{n,m}$ 。

当 $k_1\le n\le N-k_2-1$ 时，在区间 $\left[t_n,t_{n+1}\right]$ 上的配置多项式上为：

$u_h\left(t_n+sh\right)={\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1,k_2}\left(s\right)y_{n+k}}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\psi }_j\left(s\right)U_{n,j}},s\in \left[0,1\right].$

此时所使用的配置点分别为 $y_{n-k_1},y_{n-k_1+1},\cdots ,y_{n+k_2+1}$ 和 $U_{n,1},U_{n,2},\cdots ,U_{n,m}$ 。

当 $N-k_2\le n\le N-1$ 时，在区间 $\left[t_n,t_{n+1}\right]$ 上的配置多项式上为：

$u_h\left(t_n+sh\right)={\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1,k_2}\left(s\right)y_{N-k_2-1+k}}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\psi }_j\left(s\right)U_{n,j}},s\in \left[0,1\right].$

此时向后取到的配置点个数不足，为了保证取到相应个数的配置点，我们向前借用前面区间端点为配置点。这样得到的配置点分别为 $y_{N-k_1-k_2-1},y_{N-k_1-k_2},\cdots ,y_N$ 和 $U_{n,1},U_{n,2},\cdots ,U_{n,m}$ 。这样可以满足在每个区间 $\left[t_n,t_{n+1}\right]$ 的配置多项式都可以取到了 $k_1+k_2+m+2$ 个配置点。其中， $y_n:=u_h\left(t_n\right)$ 表示 $u\left(t_n\right)$ 的近似值， $U_{n,j}:=u_h\left(t_{n,j}\right)$ 表示 $u\left(t_{n，j}\right)$ 的近似值。令上述配置多项式在 $t_{n,i}$ 点精确满足方程(2)，且 $y_{n+1}=u_h\left(t_{n+1}\right)$ ，有

$\{\begin{array}{l}{U}_{n,i}=F_{n,i}+{\Phi }_{n,i},\\ y_{n+1}={\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1,k_2}}\left(1\right)y_{n+k}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\psi }_j\left(1\right)U_{n,j}.}\end{array}$

方程(1)的近似解 $u_h$ 满足下列配置方程

$f\left(t_n\right)={\displaystyle {\int }_0^{t_n}K}\left(t,s\right)u_h\left(s\right)\text{d}s,t_n\in I_h.$

将配置多项式分别带入配置方程，可以将第二类Volterra积分方程转化为一个线性方程组，求解这个线性方程组，就可以得到网格上的配置解。这个过程可以大大减少计算量和计算时间，提高数值求解的效率和准确性。

为了将上述线性系统表现为矩阵形式，定义符号 ${\alpha }_{a,b,k}^c,{\beta }_{a,b,k},{\gamma }_{a,b,j}^c,{\eta }_{a,b,j}$ ，分别有如下表达式：

${\alpha }_{a,b,k}^c={\displaystyle {\int }_0^{1}K}\left(t_{a,c},t_b+sh\right){\varphi }_k^{k_1,k_2}\left(s\right)\text{d}s,$

${\beta }_{a,b,k}={\displaystyle {\int }_0^{c_b}K}\left(t_{a,b},t_a+sh\right){\varphi }_k^{k_1,k_2}\left(s\right)\text{d}s,$

${\gamma }_{a,b,j}^c={\displaystyle {\int }_0^{1}K}\left(t_{a,c},t_b+sh\right){\psi }_j\left(s\right)\text{d}s,$

${\eta }_{a,b,j}={\displaystyle {\int }_0^{c_b}K}\left(t_{a,b},t_a+sh\right){\psi }_j\left(s\right)\text{d}s.$

我们有

$\begin{array}{c}{A}_1=\left(a_{i,j}\right)\left(i=a\left(m+1\right)+c,j=b\left(m+1\right)+1\right)\\ =\{\begin{array}{ll}{\alpha }_{a,d,b-k_1}^c,\hfill & a=1,\cdots ,k_1;b=0,\cdots ,k_1+k_2+1;c=0,\cdots ,m;d=0,\cdots ,a-1,\hfill \\ {\alpha }_{a,d,b-d}^c,\hfill & a=k_1+1,\cdots ,N-k_2;b=d-k_1,\cdots ,d+k_2+1;c=0,\cdots ,m;d=0,\cdots ,a-1,\hfill \\ {\alpha }_{a,d,b-N+k_2}^c,\hfill & a=N-k_2+1,\cdots ,N;b=N-k_1-k_2-1,\cdots ,N;c=0,\cdots ,m;d=0,\cdots ,a-1,\hfill \\ 0,\hfill & \text{others}.\hfill \end{array}\end{array}$

$\begin{array}{c}{A}_2=\left(b_{i,j}\right)\left(i=a\left(m+1\right)+c,j=b\left(m+1\right)+1\right)\\ =\{\begin{array}{ll}{\beta }_{a,b,b-k_1},\hfill & a=0,\cdots ,k_1-1;b=0,\cdots ,k_1+k_2+1;c=0,\cdots ,m,\hfill \\ {\beta }_{a,b,b-a},\hfill & a=k_1,\cdots ,N-k_2-1;b=a-k_1,\cdots ,a+k_2+1;c=0,\cdots ,m,\hfill \\ {\beta }_{a,b,b-N+k_2},\hfill & a=N-k_2,\cdots ,N-1;b=N-k_1-k_2-1,\cdots ,N;c=0,\cdots ,m,\hfill \\ 0,\hfill & \text{others}.\hfill \end{array}\end{array}$

$\begin{array}{c}{A}_3=\left(c_{i,j}\right)\left(i=a\left(m+1\right)+c,j=b\left(m+1\right)+1+d\right)\\ =\{\begin{array}{ll}{\gamma }_{a,b,j}^c,\hfill & a=1,\cdots ,N;b=0,\cdots ,N-1;c=0,\cdots ,m;d=1,\cdots ,m,\hfill \\ 0,\hfill & \text{others}.\hfill \end{array}\end{array}$

$\begin{array}{c}{A}_4=\left(d_{i,j}\right)\left(i=a\left(m+1\right)+b,j=a\left(m+1\right)+1+c\right)\\ =\{\begin{array}{ll}{\eta }_{a,b,c},\hfill & a=0,1,\cdots ,N-1;b=1,\cdots ,m;c=1,\cdots ,m,\hfill \\ 0,\hfill & \text{others}.\hfill \end{array}\end{array}$

其中 $a\_\left\{i,j\right\},b\_\left\{i,j\right\},c\_\left\{i,j\right\}$ 和 $d\_\left\{i,j\right\}$ 分别表示矩阵 $A_1,A_2,A_3$ 和 $A_4$ 第i行第j列元素。且 $1\le i\le N\left(m+1\right)$ ， $1\le j\le N\left(m+1\right)+1$ 。矩阵 $A_1,A_2,A_3$ 和 $A_4$ 均为 $N\left(m+1\right)$ 行， $N\left(m+1\right)+1$ 列的矩阵。配置方程矩阵形式：

$\left(I-hA\left(1:N\left(m+1\right),2:N\left(m+1\right)+1\right)\right)Y=F+h{y}_{0}A\left(1:N\left(m+1\right),1:1\right),$

这里， $I$ 表示单位矩阵， $A=A_1+A_2+A_3+A_4$ 。

$Y={\left(Y_1,Y_2,\cdots ,Y_N\right)}^{\text{T}},Y_n={\left(U_{n-1,0},\cdots ,U_{n-1,m},y_n\right)}^{\text{T}},n=1,\cdots ,N$

$F={\left(F_1,F_2,\cdots ,F_N\right)}^{\text{T}},F_n={\left(f\left(t_{n-1,1}\right),\cdots ,f\left(t_{n-1,m}\right),f\left(t_{n,0}\right)\right)}^{\text{T}},n=1,\cdots ,N$

3. 数值实验

在本节中，我们通过一个的数值实验来验证第二类Volterra积分方程广义多步配置法的有效性。本文所有数值实验都是在MATLAB中实现的，在分析过程中，我们用绝对误差的无穷大范数作为衡量标准，以此描述广义多步配置方法的收敛速度。观察在选取不同配置点 $k_1,k_2$ 和m时，随着N的增加，数值方法绝对误差的变化情况以及其相应的收敛阶。在这里，误差error为配置解与精确解之间绝对误差的无穷范数，收敛阶由

${\log }_2\frac{{\text{error}}_{N_{p_1}}}{{\text{error}}_{N_{p_2}}}/{\log }_2\frac{N_{p_1}}{N_{p_2}}.$

来计算。

例：考虑如下第二类Volterra积分方程

${\displaystyle {\int }_2^{t}2\cos \left(t-s\right)u\left(s\right)\text{d}s}=u\left(t\right),t\in \left[0,5\right]$

该方程的准确解为

$u\left(t\right)={\left(1+t\right)}^2{\text{e}}^t.$

我们采用了等距划分的方式对给定区间 $\left[0,5\right]$ 进行了细致的划分，通过广义多步配置法求解该方程，分别选取不同配置点 $k_1,k_2$ 和m。为了直观地展现求解结果的精确性和收敛速度，我们将绝对误差和收敛阶分别在表1和表2中列出。通过分析表格中的数据可以得到，广义多步配置法具有较快的收敛速度。

Table 1. The absolute error and convergence order of S G M C k 1 - m - k 2 ( k 1 < k 2 ) solving the above equation when m = 1

表1. 当 $m=1$ 时 $SGM{C}_{k_1-m-k_2}\left(k_1<k_2\right)$ 求解上述方程产生的绝对误差和收敛阶

在表1中，我们列出了几种 $m=1$ ，即在区间 $\left[t_n,t_{n+1}\right]$ 中间取1个点时的误差和收敛阶。可以看出，当在区间 $\left[t_n,t_{n+1}\right]$ 中间取1个配置点，观察误差的变化情况和收敛阶的观察，我们可以发现：当向后取的配置点个数 $k_2$ 大于向前取的配置点个数 $k_1$ 时，即向后取点个数大于向前取点个数 $\left(k_1<k_2\right)$ 的情况下，收敛阶会随着向后取点个数的增加而增大，分别可以达到5、6、7阶。这一现象表明，广义多步配置法在解决此类问题时具有较高的精度和收敛效率，具有重要的应用价值。

Table 2. The absolute error and convergence order of S G M C k 1 - m - k 2 ( k 1 ≥ k 2 ) solving the above equation when m = 1

表2. 当 $m=1$ 时 $SGM{C}_{k_1-m-k_2}\left(k_1\ge k_2\right)$ 求解上述方程产生的绝对误差和收敛阶

在表2中，我们列出了 $m=1$ 时，向前取点个数大于等于向后取点个数，即 $k_1\ge k_2$ 情况下收敛阶的变化情况，观察其误差的变化情况以及其相应的收敛阶。发现随着 $k_1$ 的增大，即向前取点的个数增加时，收敛阶呈现出明显的上升趋势，分别可以达到6、7、8阶。这些数据表明，随着向前取点个数的增加，广义多步配置法的收敛速度越来越快。

Table 3. The absolute error and convergence order of S G M C k 1 - m - k 2 solving the above equation when m = 2

表3. 当 $m=2$ 时 $SGM{C}_{k_1-m-k_2}$ 求解上述方程产生的绝对误差和收敛阶

在表3中，我们列出了 $m=2$ 时的几种情况，即在小区间中取两个配置点，分别改变向前向后取点个数，观察收敛阶变化，我们发现，取点个数越多，收敛阶越高。这些数据表明，随着取点个数的增加，广义多步配置法的收敛速度加快。

通过对表1~3中的数据进行分析可以得到：广义多步配置法随着N取值的增加，广义多步配置法的误差表现呈现出逐渐减小的趋势，且随着配置点个数 $k_1+k_2+m$ 的增加，该方法的收敛阶也在稳步上升，且其收敛阶可以达到 $k_1+k_2+m+2$ 阶，说明用广义多步配置法求解第二类Volterra积分方程是有效的。

4. 总结

本文在经典多步配置法的基础上，引入了边值方法的思想，创新性地提出了针对第二类Volterra积分方程的广义多步配置法。首先，我们构建了一种新的广义多步配置法，将原方程进行精细的离散化处理，将其转化为一组线性方程组。通过求解该线性系统，我们可以得到配置解。接下来，为了验证此类方法的可行性和有效性，我们进行了详细的数值实验。

实验结果表明，广义多步配置法在求解第二类Volterra积分方程时具有显著的优势，相对于经典的多步配置方法，可以达到更高的收敛阶，突显其在求解第二类Volterra积分方程中的重要性。