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

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

摘要: Volterra积分方程广泛应用在许多科学研究领域，例如热传导模型、声学散射问题、人口预测模型等。针对第一类Volterra积分方程的数值解，在经典多步配置法的基础上利用边值方法的思想，研究其广义多步配置方法。利用Lagrange插值公式，选取不同的节点作为插值节点，将原方程离散成为一个线性方程组。通过实验验证了该方法求解第一类Volterra积分方程的有效性，并且该方法可以达到较高的收敛阶。

Abstract: The Volterra integral equation is widely used in many scientific research fields, such as heat transfer models, acoustic scattering problems, population prediction models and so on. Aiming at the numerical solution of the first type of Volterra integral equation, the idea of edge value method is used to study its generalized multistep collocation method based on the classical multistep collocation method. Using the Lagrange interpolation formula, different nodes are selected as interpolation nodes to discretize the original equation into a linear equation system. The effectiveness of the method in solving the first type of Volterra integral equation is verified by experiments, and the method can reach a higher convergence order.

1. 引言

Volterra积分方程在物理学、生物学、化学、经济学等多个领域都有着广泛应用。常见的热传导模型、流行病扩散、声学散射问题等都可以用Volterra积分方程来模拟 [1] [2] [3] [4] 。例如种群增长模型以及疾病传播模型可以由如下形式的Volterra积分方程描述：

$u\left(t\right)=f\left(t\right){\displaystyle {\int }_0^{t}P\left(t-s\right)G\left(u\left(s\right)\right)\text{d}s}.$

对于一般的Volterra积分方程，真实解难以通过解析方法得到，所以依靠数值方法求解Volterra积分方程的近似解成为一个重要的研究方向。近年来，许多学者对数值求解Volterra积分方程展开了深入的研究，提出了许多有效的数值方法，如线性多步法、Runge-Kutta方法、谱方法、配置法等。

1987年，McAlevey基于矩形求积法则和中点求积法则求解第一类Volterra积分方程，并讨论了其误差的渐进展开式 [5] 。2004年，Brunner的专著介绍了第一类、第二类Volterra积分方程在不连续多项式空间的配置方法及其收敛分析 [6] 。2009年，Brunner和Davies等给出了一类不连续Galerkin方法，并证明了该方法的收敛性和高振荡超收敛性 [7] 。2011年，Wang和Xiang通过分析含Bessel函数的振荡积分的渐进性质，推导了Volterra积分方程精确解的渐进展开式，并通过Filon型方法来求解其近似解 [8] 。同年，Chen和Zhang构造了Volterra积分方程的边值方法 [9] 。随后，2012年，Chen和Zhang对这类算法进行了全面的研究 [10] 。2013年，Xiang和Brunner提出了高效的Filon配置方法用于求解Volterra积分方程 [11] 。2014年，Xiang基于Laplace变换和Laplace逆变换推导了Volterra积分方程的显示表达式，并利用Clenshaw-Curtis-Filon型方法来得到其近似解 [12] 。2015年，Fazeli和Hojjati介绍了Volterra积分微分方程的超隐式多步配置方法 [13] 。2018年，Xiang和Li等讨论了求解滞后势能积分方程的有效数值方法，通过连续Fourier变换将滞后势能积分方程转化为Volterra积分方程，然后利用Clenshaw-Curtis-type方法来求解该方程得到其近似解 [14] 。2019年，Li和Xiang等研究了带有高振荡核的第一类Volterra积分方程，利用Laplace变换和Laplace逆变换推导了Volterra积分方程的另一显示表达式，并使用Clenshaw-Curtis-Filonh-type和Clenshaw-Curtis-type等方法求解方程得到其近似解 [15] 。

用配置法求解Volterra积分方程不仅可以获取较好的数值精度，同时计算量相对较小。多步配置法由于收敛性好，易于构造，广泛受到国内外学者的重视。2009年，Conte和Paternoster利用单步方法得到的近似值构造了多步配置方法。该方法能够在不增加配置点的情况下提高多步配置方法的收敛阶，进一步分析了该方法的收敛性及其线性稳定性 [16] 。2012年，Liang和Brunner研究了第一类Volterra积分方程的配置方法，并对其配置解的局部超收敛性进行了分析 [17] 。2017年，Ma和Xiang基于特殊的多步配置方法，利用未计算的近似值，提出求解Volterra积分方程的配置边值方法，并分析了其线性稳定性 [18] 。2018年，Zhang和Liang将多步配置法应用于第一类Volterra积分方程，证明了多步配置解的存在唯一性，并分析了该方法的收敛条件，给出了相应的收敛阶 [19] 。2019年，Zhao和Long等人利用光滑变换将弱奇异Volterra积分方程转化为具有更好正则性的弱奇异Volterra积分方程，变换后的积分方程的解充分光滑。然后利用已计算的近似值和当前以及下一个子区间的配置点构造了超隐式多步配置法，并分析了其收敛性和稳定性 [20] 。2021年，Liu和Ma研究了第一类Volterra积分方程的一类块配置边值方法。通过研究配置方程的特殊结构，讨论了其可解性，并给出了配置解存在的充分条件，通过数值实验验证了方法的有效性 [21] 。2022年，Patil和Shinde等利用Anuj变换来求解第一类线性Volterra积分方程。证明了Anuj变换对于求解这类方程是有效的 [22] 。2022年，Zhao和Fan等针对带有高振荡核的线性Volterra积分方程提出配置方法，采用了Filon型方法对配置方程中的振荡积分进行离散，并研究了该方法的收敛性 [23] 。

当使用配置法求解Volterra积分方程时，可以得到数值解对真实解的整体逼近，且具有较好的绝对稳定域。本文对第一类Volterra积分方程提出了广义多步配置方法，用于解决大规模的复杂问题。

2. 广义多步配置法的构造

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

$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(t)={\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.$

具体地，当 $n=1,\cdots ,k_1$ ，有

$\begin{array}{l}f\left(t_{n,i}\right)=F_{n,i}+{\Phi }_{n,i}\\ =h{\displaystyle \underset{l_1=0}{\overset{n-1}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_1}+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{n,j}\right)\text{d}s}}\\ +h{\displaystyle {\int }_0^{c_i}k\left(t_{n,i},t_n+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{n,j}\right)\text{d}s},\end{array}$

当 $n=k_1+1,\cdots ,N-k_2$ ，有

$\begin{array}{l}f\left(t_{n,i}\right)=F_{n,i}+{\Phi }_{n,i}\\ =h{\displaystyle \underset{l_1=0}{\overset{k_1-1}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_1}+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{k_1,j}\right)\text{d}s}}\\ +h{\displaystyle \underset{l_2=k_1+1}{\overset{n}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_2-1}+sh\right)\left({\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1{k}_2}}\left(s\right)y_{l_2-1+k}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\phi }_j}\left(s\right)U_{l_2-1,j}\right)\text{d}s}}\\ +h{\displaystyle {\int }_0^{c_i}k\left(t_{n,i},t_n+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{n,j}\right)\text{d}s},\end{array}$

当 $n=N-k_2+1,\cdots ,N$ ，有

$\begin{array}{l}f\left(t_{n,i}\right)=F_{n,i}+{\Phi }_{n,i}\\ =h{\displaystyle \underset{l_1=0}{\overset{k_1-1}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_1}+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{k_1,j}\right)\text{d}s}}\\ +h{\displaystyle \underset{l_2=k_1+1}{\overset{N-k_2}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_2-1}+sh\right)\left({\displaystyle \underset{k=-k_1}{\overset{k_2+1}{\sum }}{\varphi }_k^{k_1{k}_2}}\left(s\right)y_{l_2-1+k}+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\phi }_j}\left(s\right)U_{l_2-1,j}\right)\text{d}s}}\\ +h{\displaystyle \underset{l_3=N-k_2+1}{\overset{N-1}{\sum }}{\displaystyle {\int }_0^{1}k\left(t_{n,i},t_{l_3}+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{l_3,j}\right)\text{d}s}}\\ +h{\displaystyle {\int }_0^{c_i}k\left(t_{n,i},t_n+sh\right)\left({\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 }}{\phi }_j}\left(s\right)U_{n,j}\right)\text{d}s.}\end{array}$

通过上述过程，就实现了将第一类Volterra积分方程离散为一个线性系统，求解该线性系统，就可以得到网格上的配置解。

3. 数值实验

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

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

来计算。

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

${\displaystyle {\int }_0^{t}2\cos \left(t-s\right)u\left(s\right)\text{d}s=\left({\left(1+t\right)}^2-1\right)e^t,}t\in \left[0,4\right],$

该方程的准确解为

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

对区间 $\left[0,4\right]$ 进行等距划分，通过广义多步配置法求解该方程，分别选取不同配置点 $k_1,k_2$ 和m，绝对误差和收敛阶在表1和表2中列出。通过分析表格中的数据可以得到，广义多步配置法具有较快的收敛速度。

Table 1. Absolute errors and convergence orders of F G M C k 1 - 0 - k 2 for Example 1

表1. $FGM{C}_{k_1\text{-}0\text{-}{k}_2}$ 求解例1时产生的绝对误差和收敛阶

Table 2. Absolute errors and convergence orders of F G M C k 1 - m - k 2 for Example 1

表2. $FGM{C}_{k_1-m-k_2}$ 求解例1时产生的绝对误差和收敛阶

在表1中，我们列出了几种 $m=0$ ，即在区间 $\left[t_n,t_{n+1}\right]$ 中间不取点时的误差和收敛阶。可以看出，当在区间 $\left[t_n,t_{n+1}\right]$ 中间不取配置点，且向后取的配置点个数 $k_2$ 大于等于向前取的配置点个数 $k_1$ 时， $FGM{C}_{k_1-0-k_2}$ 的收敛阶可以达到 $k_1+k_2+2$ 阶。

在表2中，我们列出了 $m\ne 0$ 时的几种情况，观察其误差的变化情况以及其相应的收敛阶。发现当 $m=1$ 时，向前取1个配置点( $k_1=1$ )，向后取3个配置点( $k_2=3$ )时，收敛阶可以达到7阶。当 $m=2$ 时，向前取1个配置点( $k_1=1$ )，向后取2个配置点( $k_2=2$ )时，收敛阶可以达到7阶，向前取2个配置点( $k_1=2$ )，向后取2个配置点( $k_2=2$ )时，收敛阶可以达到8阶。

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

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

该方程的准确解为

$u\left(t\right)=2\cos t-1.$

同样使用广义多步配置法求解此方程，绝对误差和收敛阶的实验结果在表3和表4中列出。实验结果验证了该方法可以达到高收敛阶，在表3中可以看出当 $m=0$ 且 $k_1\le k_2$ 时，收敛阶可以达到 $k_1+k_2+2$ 阶。由表4得到，当 $m=1$ 时，取 $k_1=1,k_2=3$ 收敛阶可以达到7阶。当 $m=2$ 时，取 $k_1=1,k_2=2$ 和 $k_1=2,k_2=2$ 收敛阶分别可以达到7阶和8阶。

Table 3. Absolute errors and convergence orders of F G M C k 1 - 0 - k 2 for Example 2

表3. $FGM{C}_{k_1-0-k_2}$ 求解例2时产生的绝对误差和收敛阶

Table 4. Absolute errors and convergence orders of F G M C k 1 - m - k 2 for Example 2

表4. $FGM{C}_{k_1-m-k_2}$ 求解例2时产生的绝对误差和收敛阶

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

4. 总结

本文在经典多步配置法的基础上利用边值方法的思想，提出了第一类Volterra积分方程的广义多步配置法 $FGM{C}_{k_1-m-k_2}$ 。首先构造广义多步配置法方法，将原方程离散成一个线性系统，通过求解线性系统得到配置解。然后通过数值实验表明，此类方法是有效的并且相对于经典的多步配置方法可以达到较高的收敛阶。

基金项目

国家自然科学基金项目(11901133)。

NOTES

^*通讯作者。