指标-3型积分代数方程的配置边值方法

Collocation Boundary Value Method for Index-3 Integral Algebraic Equations

摘要: 针对指标-3型积分代数方程的数值解，研究其配置边值方法，基于插值多项式，利用未计算的近似值，通过将原方程进行离散化构造了指标-3型积分代数方程的配置边值方法，并分析了该方法的可解性和收敛性，证明了利用该方法求解指标-3型积分代数方程可达到较高收敛阶，最后通过数值实验验证了方法的有效性。

Abstract: Regarding the numerical solution of the index-3 integral algebraic equation, the collocation boundary value method was investigated. Based on the interpolation polynomial and the utilization of uncomputed approximate values, the collocation boundary value method for the index-3 integral algebraic equation was constructed by discretizing the original equation. The solvability and convergence of this method were analyzed. It was demonstrated that the application of this method in solving the index-3 integral algebraic equation can achieve a relatively high convergence order. Finally, the validity of the method was verified through numerical experiments.

1. 引言

积分代数方程是第一类和第二类Volterra积分方程的混合系统，该方程有许多实际应用，比如粘弹性材料领域[1] [2]、热传导过程中存储器核的识别[3]以及在一个小细胞内的化学反应的进化[4]。该方程的一般形式如下

$A\left(t\right)X\left(t\right)=G\left(t\right)+{\displaystyle {\int }_0^{t}K\left(t,s,X\left(s\right)\right)\text{d}s},t\in I=\left[0,T\right],$

其中$X,G:I\to {ℝ}^d,K:I\times I\times {ℝ}^d\to {ℝ}^d$ 在定义域上是连续的，并且$A\left(t\right)\in L\left({ℝ}^d\right)$ 是奇异矩阵且$rank\left(A\right)>0$ 。有许多问题，都可以用微分代数方程以以下形式表示：

$A\left(t\right)X^{\prime }\left(t\right)=F\left(t,X\left(t\right)\right),$

其中$F:I\times {ℝ}^d\to {ℝ}^d$ 。通过观察可知，微分代数方程可以通过积分转换为积分代数方程，进一步，指标-3型微分代数方程就可以转化为指数-3型积分代数方程[5]。由于积分代数方程的解析解很难直接求出，因此寻找合适的数值方法求解积分代数方程就变得十分重要。

考虑第一类和第二类Volterra积分方程的混合系统：

$\{\begin{array}{l}x\left(t\right)=f\left(t\right)+{\displaystyle {\int }_0^t{K}_{11}\left(t,s\right)x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{12}\left(t,s\right)y\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{13}\left(t,s\right)z\left(s\right)\text{d}s},\\ y\left(t\right)=g\left(t\right)+{\displaystyle {\int }_0^t{K}_{21}\left(t,s\right)x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{22}\left(t,s\right)y\left(s\right)\text{d}s},\\ 0=h_1\left(t\right)+{\displaystyle {\int }_0^t{K}_{32}\left(t,s\right)y\left(s\right)\text{d}s},t\in I=\left[0,T\right],\end{array}$ (1)

其中

$x,f:I\to {ℝ}^{d_1},y,g:I\to {ℝ}^{d_2},$ $z,h_1:I\to {ℝ}^{d_3}.$

矩阵核

$K_{kk}\left(\cdot ,\cdot \right)\in L\left({ℝ}^{d_k}\right),K_{12}\left(\cdot ,\cdot \right)\in L\left({ℝ}^{d_2},{ℝ}^{d_1}\right),K_{13}\left(\cdot ,\cdot \right)\in L\left({ℝ}^{d_3},{ℝ}^{d_1}\right),$

$K_{21}\left(\cdot ,\cdot \right)\in L\left({ℝ}^{d_1},{ℝ}^{d_2}\right),K_{32}\left(\cdot ,\cdot \right)\in L\left({ℝ}^{d_2},{ℝ}^{d_3}\right)$

是连续的，并且$L\left(\cdot ,\cdot \right)$ 是线性变换空间。

我们将系统(1)称为指标-3型的半显式积分代数方程，因为它需要三次微分才能将其简化为一个正则的积分方程组。然而，从数值的角度来看，将它简化为指标-2型的积分代数方程是不实际的。

迄今为止，已有诸多学者研究了积分代数方程的数值解[5]-[22]。2000年，Kauthen [6]研究了指标-1型积分代数方程的多项式样条配置方法，并给出了全局收敛性和局部超收敛性分析。2020年，Zhang等[7]引入了多步配置方法来求解指标-1型积分代数方程，证明了其多步配置解的存在唯一性，研究了摄动多步配置方法的收敛性，并扩展了包括对无摄动项的多步配置方法的分析。2013年，Bulatov和Budnikova [8]测试了指标-1型积分代数方程的一般多步法的稳定性。2018年，Farahani和Hadizadeh [9]用一种直接正则化方法，对指标-1型积分代数方程的数值解构造了一种基于Lavrentiev的正则化迭代方法的数值算法，并对该方法的收敛性进行了分析。2013年，Liang和Brunner [10]基于Volterra积分算子的v-平滑性质的概念，将给定的积分代数方程系统解耦成正则的第二类Volterra积分方程的固有系统和一个第一类Volterra积分方程系统，然后利用该解耦方法，得到了指标-1型积分代数方程分段多项式配置解的最优收敛性。2016年，Liang和Brunner [11]研究了可处理性指标$\mu =2$ 和$\left(v+1\right)$ -平滑$\left(v\ge 1\right)$ 的积分代数方程系统的解耦，然后证明了配置方法在给定的该系统中的应用与在解耦系统中的应用是等价的，并对半显式指标-2型积分代数方程的配置解进行了收敛性分析。2011年，Hadizadeh等[12]针对指标-2型积分代数方程，提出了一种包含矩阵–向量乘法表示的Jacobi配置方法，对加权$L^2$ 范数中的误差界进行了严格的分析，从理论上证明了当核和原函数足够光滑时的谱收敛速度。2012年，Ghoreishi等[13]研究了指标-2型积分代数方程的多项式样条配置方法，给出了保证该方法收敛的充要条件，还分析了配置参数$c_m$ 的两种不相交情况下的收敛速度。2015年，Pishbin [5]研究了指标-3型积分代数方程的分段多项式配置方法，分析了该配置方法的全局收敛性，并建立了收敛结果的最优阶数。

另一方面，配置边值方法是求解Volterra积分方程的一种有效方法。2017年，Ma和Xiang [23]研究了Volterra积分方程的边值方法，利用特殊的多步配置方法设计了高阶数值格式，该方法依赖于接下来几步中解的数值逼近，稳定性分析表明，这些方法具有广泛的绝对稳定区域。2019年，Ma和Liu [24]针对一类弱奇异Volterra积分方程的数值解，提出了有效的分数配置边值方法，并推导出了局部收敛估计。2022年，Liu和Ma [25]分析了第二类自卷积Volterra积分方程的配置边值解。

本文利用配置边值方法来数值求解指标-3型积分代数方程(1)，其主要工作安排如下，在第一节中，利用Lagrange基函数来构造指标-3型积分代数方程的配置边值方法。在第二节中，利用数值例子来验证该方法的有效性。最后进行总结。

2. 配置边值方法

Ma和Xiang [23]在2017年提出了线性Volterra积分方程的配值边值方法，他们利用了特殊的多步配置方法提出的这类边值方法，可以看作是多步配置方法的扩展。在这里，我们将这类边值方法应用于指标-3型积分代数方程。

首先对区间$\left[0,T\right]$ 作等距网格划分

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

定义$X_h=I_h\backslash \left\{0\right\}$ ，在区间$\left[t_n,t_{n+1}\right]$ 上，选取$t_n$ 后的k个节点为配置点，定义基函数为

${\varphi }_j^k\left(s\right)={\displaystyle \prod _{i=0,i\ne j}^{k+1}\frac{s-i}{j-i}},$

则定义区间$\left[t_n,t_{n+1}\right]\left(n=0,1,\cdots ,N-k-1\right)$ 上的配置多项式$x_h,y_h,z_h$ 分别为

$\{\begin{array}{l}{x}_h\left(t_n+sh\right)={\displaystyle \sum _{i=0}^{k+1}{\varphi }_i^k\left(s\right)x_h}\left(t_{n+i}\right),\\ y_h\left(t_n+sh\right)={\displaystyle \sum _{i=0}^{k+1}{\varphi }_i^k\left(s\right)y_h}\left(t_{n+i}\right),\\ z_h\left(t_n+sh\right)={\displaystyle \sum _{i=0}^{k+1}{\varphi }_i^k\left(s\right)z_h}\left(t_{n+i}\right),\end{array}$ (2)

在区间$\left[t_{N-k},t_N\right]$ 上选取该区间的全部插值节点(共$k+1$ 个)来构造配置多项式，定义基函数为

${\varphi }_j^{k-1}\left(s\right)={\displaystyle \prod _{i=0,i\ne j}^k\frac{s-i}{j-i}},$

定义该区间上的配置多项式$x_h,y_h,z_h$ 分别为

$\{\begin{array}{l}{x}_h\left(t_{N-k}+sh\right)={\displaystyle \sum _{i=0}^k{\varphi }_i^{k-1}\left(s\right)x_h}\left(t_{N-k+i}\right),\\ y_h\left(t_{N-k}+sh\right)={\displaystyle \sum _{i=0}^k{\varphi }_i^{k-1}\left(s\right)y_h}\left(t_{N-k+i}\right),\\ z_h\left(t_{N-k}+sh\right)={\displaystyle \sum _{i=0}^k{\varphi }_i^{k-1}\left(s\right)z_h}\left(t_{N-k+i}\right),\end{array}$ (3)

其中，$x_h\left(t_n\right),y_h\left(t_n\right),z_h\left(t_n\right)$ 分别表示$x\left(t_n\right),y\left(t_n\right),z\left(t_n\right)$ 的近似值。因此，就可以得到方程组(1)的k-步配置边值方法(k-CBVM)的配置方程：

$\{\begin{array}{l}{x}_h\left(t\right)=f\left(t\right)+{\displaystyle {\int }_0^t{K}_{11}\left(t,s\right)x_h\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{12}\left(t,s\right)y_h\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{13}\left(t,s\right)z_h\left(s\right)\text{d}s},\\ y_h\left(t\right)=g\left(t\right)+{\displaystyle {\int }_0^t{K}_{21}\left(t,s\right)x_h\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{22}\left(t,s\right)y_h\left(s\right)\text{d}s},\\ 0=h_1\left(t\right)+{\displaystyle {\int }_0^t{K}_{32}\left(t,s\right)y_h\left(s\right)\text{d}s},t\in X_h.\end{array}$

通过求解上述线性系统，就可以得到方程组(1)的配置边值解。

特别地，对于1-步配置边值方法(1-CBVM)，在子区间$\left[t_n,t_{n+1}\right]\left(n<N-1\right)$ 上的配置解$x_h,y_h,z_h$ 分别为

$\{\begin{array}{l}{x}_h\left(t_n+sh\right)=x_h\left(t_n\right){\varphi }_0^1\left(s\right)+x_h\left(t_{n+1}\right){\varphi }_1^1\left(s\right)+x_h\left(t_{n+2}\right){\varphi }_2^1\left(s\right),\\ y_h\left(t_n+sh\right)=y_h\left(t_n\right){\varphi }_0^1\left(s\right)+y_h\left(t_{n+1}\right){\varphi }_1^1\left(s\right)+y_h\left(t_{n+2}\right){\varphi }_2^1\left(s\right),\\ z_h\left(t_n+sh\right)=z_h\left(t_n\right){\varphi }_0^1\left(s\right)+z_h\left(t_{n+1}\right){\varphi }_1^1\left(s\right)+z_h\left(t_{n+2}\right){\varphi }_2^1\left(s\right),\end{array}$

在子区间$\left[t_{N-1},t_N\right]$ 上的配置解$x_h,y_h,z_h$ 分别为

$\{\begin{array}{l}{x}_h\left(t_{N-1}+sh\right)=x_h\left(t_{N-1}\right){\varphi }_0^0\left(s\right)+x_h\left(t_N\right){\varphi }_1^0\left(s\right),\\ y_h\left(t_{N-1}+sh\right)=y_h\left(t_{N-1}\right){\varphi }_0^0\left(s\right)+y_h\left(t_N\right){\varphi }_1^0\left(s\right),\\ z_h\left(t_{N-1}+sh\right)=z_h\left(t_{N-1}\right){\varphi }_0^0\left(s\right)+z_h\left(t_N\right){\varphi }_1^0\left(s\right).\end{array}$

定义

${\alpha }_{pq,i,j}^{\tau }={\displaystyle {\int }_0^1{K}_{pq}\left(ih,jh+sh\right){\varphi }_{\tau }^1\left(s\right)\text{d}s}$

和

${\beta }_{pq,i,j}^{\upsilon }={\displaystyle {\int }_0^1{K}_{pq}\left(ih,jh+sh\right){\varphi }_{\upsilon }^0\left(s\right)\text{d}s}$

其中$p,q=1,2,3$ ，$i=1,2,\cdots ,N$ ，$j=0,1,\cdots ,N-1$ ，$\tau =0,1,2$ ，$\upsilon =0,1$ ，则可以构造$N\times N$ 阶矩阵：

$A_{pq}:=L_{pq}+M_{pq}+R_{pq},$

其中

$L_{pq}=\left(\begin{array}{cccccccc}0& 0& 0& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,2,1}^0& 0& 0& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,3,1}^0& {\alpha }_{pq,3,2}^0& 0& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,4,1}^0& {\alpha }_{pq,4,2}^0& {\alpha }_{pq,4,3}^0& 0& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ {\alpha }_{pq,N-1,1}^0& {\alpha }_{pq,N-1,2}^0& {\alpha }_{pq,N-1,3}^0& {\alpha }_{pq,N-1,4}^0& \cdots & {\alpha }_{pq,N-1,N-2}^0& 0& 0\\ {\alpha }_{pq,N,1}^0& {\alpha }_{pq,N,2}^0& {\alpha }_{pq,N,3}^0& {\alpha }_{pq,N,4}^0& \cdots & {\alpha }_{pq,N,N-2}^0& {\beta }_{pq,N,N-1}^0& {\beta }_{pq,N,N-1}^1\end{array}\right),$

$M_{pq}=\left(\begin{array}{cccccccc}{\alpha }_{pq,1,0}^1& 0& 0& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,2,0}^1& {\alpha }_{pq,2,1}^1& 0& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,3,0}^1& {\alpha }_{pq,3,1}^1& {\alpha }_{pq,3,2}^1& 0& \cdots & 0& 0& 0\\ {\alpha }_{pq,4,0}^1& {\alpha }_{pq,4,1}^1& {\alpha }_{pq,4,2}^1& {\alpha }_{pq,4,3}^1& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ {\alpha }_{pq,N-1,0}^1& {\alpha }_{pq,N-1,1}^1& {\alpha }_{pq,N-1,2}^1& {\alpha }_{pq,N-1,3}^1& \cdots & {\alpha }_{pq,N-1,N-3}^1& {\alpha }_{pq,N-1,N-2}^1& 0\\ {\alpha }_{pq,N,0}^1& {\alpha }_{pq,N,1}^1& {\alpha }_{pq,N,2}^1& {\alpha }_{pq,N,3}^1& \cdots & {\alpha }_{pq,N,N-3}^1& {\alpha }_{pq,N,N-2}^1& 0\end{array}\right),$

$R_{pq}=\left(\begin{array}{cccccccc}0& {\alpha }_{pq,1,0}^2& 0& 0& \cdots & 0& 0& 0\\ 0& {\alpha }_{pq,2,0}^2& {\alpha }_{pq,2,1}^2& 0& \cdots & 0& 0& 0\\ 0& {\alpha }_{pq,3,0}^2& {\alpha }_{pq,3,1}^2& {\alpha }_{pq,3,2}^2& \cdots & 0& 0& 0\\ 0& {\alpha }_{pq,4,0}^2& {\alpha }_{pq,4,1}^2& {\alpha }_{pq,4,2}^2& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& {\alpha }_{pq,N-1,0}^2& {\alpha }_{pq,N-1,1}^2& {\alpha }_{pq,N-1,2}^2& \cdots & {\alpha }_{pq,N-1,N-4}^2& {\alpha }_{pq,N-1,N-3}^2& {\alpha }_{pq,N-1,N-2}^2\\ 0& {\alpha }_{pq,N,0}^2& {\alpha }_{pq,N,1}^2& {\alpha }_{pq,N,2}^2& \cdots & {\alpha }_{pq,N,N-4}^2& {\alpha }_{pq,N,N-3}^2& {\alpha }_{pq,N,N-2}^2\end{array}\right),$

因此，配置方程可以转化为矩阵形式

$\left(\begin{array}{ccc}I+h{A}_{11}& h{A}_{12}& h{A}_{13}\\ h{A}_{21}& I_N+h{A}_{11}& O\\ O& h{A}_{32}& O\end{array}\right)\left(\begin{array}{c}{x}_h\\ \begin{array}{l}{y}_h\\ z_h\end{array}\end{array}\right)=\left(\begin{array}{c}{r}_1\\ \begin{array}{l}{r}_2\\ r_3\end{array}\end{array}\right),$ (4)

这里$I$ 表示N阶单位矩阵，$O$ 表示N阶零矩阵，

$r_1=\left(\begin{array}{c}f\left(h\right)-hx\left(0\right){\alpha }_{11,1,0}^0-hy\left(0\right){\alpha }_{12,1,0}^0-hz\left(0\right){\alpha }_{13,1,0}^0\\ f\left(2h\right)-hx\left(0\right){\alpha }_{11,2,0}^0-hy\left(0\right){\alpha }_{12,2,0}^0-hz\left(0\right){\alpha }_{13,2,0}^0\\ ⋮\\ f\left(T\right)-hx\left(0\right){\alpha }_{11,N,0}^0-hy\left(0\right){\alpha }_{12,N,0}^0-hz\left(0\right){\alpha }_{13,N,0}^0\end{array}\right),$

$r_2=\left(\begin{array}{c}g\left(h\right)-hx\left(0\right){\alpha }_{21,1,0}^0-hy\left(0\right){\alpha }_{22,1,0}^0\\ g\left(2h\right)-hx\left(0\right){\alpha }_{21,2,0}^0-hy\left(0\right){\alpha }_{22,2,0}^0\\ ⋮\\ g\left(T\right)-hx\left(0\right){\alpha }_{21,N,0}^0-hy\left(0\right){\alpha }_{22,N,0}^0\end{array}\right),r_3=\left(\begin{array}{c}{h}_1\left(h\right)-hy\left(0\right){\alpha }_{32,1,0}^0\\ h_1\left(2h\right)-hy\left(0\right){\alpha }_{32,2,0}^0\\ ⋮\\ h_1\left(T\right)-hy\left(0\right){\alpha }_{32,N,0}^0\end{array}\right)$

以及

$x_h^{\text{T}}=\left[x_h\left(h\right),x_h\left(2h\right),\cdots ,x_h\left(T\right)\right],y_h^{\text{T}}=\left[y_h\left(h\right),y_h\left(2h\right),\cdots ,y_h\left(T\right)\right],$

$z_h^{\text{T}}=\left[z_h\left(h\right),z_h\left(2h\right),\cdots ,z_h\left(T\right)\right].$

3. 收敛性分析

指标-3型积分代数方程的配置解的存在唯一性由系统(4)的系数矩阵的可逆性决定，当系数矩阵的行列式不为0时，其配置解存在且唯一，即

$\left|\begin{array}{ccc}I+h{A}_{11}& h{A}_{12}& h{A}_{13}\\ h{A}_{21}& I_N+h{A}_{11}& O\\ O& h{A}_{32}& O\end{array}\right|=h^{3N}\left|A_{13}\right|\left|A_{21}\right|\left|A_{32}\right|\ne 0,$

又由于$A_{13},A_{21},A_{32}$ 之间只有核函数不相同，所以只要它们中其中一个矩阵行列式不为0，则系数矩阵行列式不为0。这里我们主要考虑矩阵$A_{32}$ ，即考虑方程组(1)中的第三个方程。

对于方程组(1)，可以将其变形为如下形式：

$\{\begin{array}{l}-f\left(t\right)=x\left(t\right)+{\displaystyle {\int }_0^t{K}_{11}\left(t,s\right)x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{12}\left(t,s\right)y\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{13}\left(t,s\right)z\left(s\right)\text{d}s},\\ -g\left(t\right)=y\left(t\right)+{\displaystyle {\int }_0^t{K}_{21}\left(t,s\right)x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{K}_{22}\left(t,s\right)y\left(s\right)\text{d}s},\\ -h_1\left(t\right)={\displaystyle {\int }_0^t{K}_{32}\left(t,s\right)y\left(s\right)\text{d}s},t\in I=[0,T],\end{array}$ (5)

下面我们首先考虑积分代数方程(5)中核函数都等于1的特殊情况。而对于一般的核函数，其存在唯一性和收敛性可类似推导。

对于方程组(5)中的第三个方程，当核函数$K_{32}\left(t,s\right)=1$ 时，它满足

$-h_1\left(t_n\right)=h{\displaystyle \sum _{j=0}^{n-1}{\displaystyle {\int }_0^1{y}_h\left(t_j+sh\right)\text{d}s}},$

对于$n=N-1,N-2,\cdots ,0$ ，计算$-h_1\left(t_{n+1}\right)-\left(-h_1\left(t_n\right)\right)$ 有

$-h_1\left(t_{n+1}\right)-\left(-h_1\left(t_n\right)\right)=h{\displaystyle {\int }_0^1{y}_h\left(t_n+sh\right)\text{d}s},$

当$n=0,1,\cdots ,N-k-2$ 时，将配置多项式(2)代入上式有

$-h_1\left(t_{n+1}\right)-\left(-h_1\left(t_n\right)\right)=h{\displaystyle \sum _{i=0}^{k+1}{\displaystyle {\int }_0^1{\varphi }_i^k\left(s\right)\text{d}s}}{y}_h\left(t_{n+i}\right),$ (6)

当$n=N-k-1,N-k,\cdots ,N-1$ 时，将配置多项式(3)代入有

$-h_1\left(t_{n+1}\right)-\left(-h_1\left(t_n\right)\right)=h{\displaystyle \sum _{i=0}^k{\displaystyle {\int }_0^{n-N+k+2}{\varphi }_i^{k-1}\left(s\right)ds}}{y}_h\left(t_{N-k+i}\right),$ (7)

为了将上述线性系统(6)和(7)表达为矩阵形式，下面先定义一些符号。令$a_m\left(t\right)={\displaystyle \sum _{i=-m+1}^{m-1}{a}_i{t}^i}$ ，$m\times m$ 阶Toeplitz矩阵$T\left[a_m\right]$ 定义为

$T\left[a_m\right]=\left(\begin{array}{cccc}{a}_0& a_{-1}& \cdots & a_{-m+1}\\ a_1& a_0& \cdots & a_{-m+2}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m-1}& a_{m-2}& \cdots & a_0\end{array}\right).$

分别定义$d_i^k$ 和${\widehat{d}}_{n,i}^k$ ：

$d_i^k:={\displaystyle {\int }_0^1{\varphi }_i^k\left(s\right)\text{d}s},{\widehat{d}}_{n,i}^k:={\displaystyle {\int }_0^{n-N+k+2}{\varphi }_i^{k-1}\left(s\right)\text{d}s}.$

Laurent多项式${\mathfrak{d}}_m^i\left(t\right)$ 定义为

${\mathfrak{d}}_m^i\left(t\right):={\displaystyle \sum _{j=-m+1}^0{d}_{j-1,i}}{t}^j,i=0,1,\cdots ,k+1.$

令矩阵$B_N^k=\left(\begin{array}{cc}{T}_N^k& r\\ 0& R_N^k\end{array}\right)$ ，其中$r$ 是一个$\left(N-k-1\right)\times \left(k+1\right)$ 阶矩阵，$T_N^k$ 是$\left(N-k-1\right)\times \left(N-k-1\right)$ 阶Toeplitz矩阵，其元素由下面的Laurent级数的系数生成：

${\text{c}}_{N-k-1}^k\left(t\right):={\displaystyle \sum _{j=-k}^1{d}_{1,1-j}^k{t}^j}$ ，其中$d_{1,1-j}^k=d_{1-j}^k={\displaystyle {\int }_0^1{\varphi }_{1-j}^k\left(s\right)\text{d}s}$ 。

$0$ 是$\left(k+1\right)\times \left(N-k-1\right)$ 阶零矩阵，

$R_N^k=\left(\begin{array}{cccc}{\widehat{d}}_{N-k-1,0}^k& {\widehat{d}}_{N-k-1,1}^k& \cdots & {\widehat{d}}_{N-k-1,k}^k\\ {\widehat{d}}_{N-k,0}^k& {\widehat{d}}_{N-k,1}^k& \cdots & {\widehat{d}}_{N-k,k}^k\\ ⋮& ⋮& \ddots & ⋮\\ {\widehat{d}}_{N-1,0}^k& {\widehat{d}}_{N-1,1}^k& \cdots & {\widehat{d}}_{N-1,k}^k\end{array}\right),$

令$H={\left(-h_1\left(t_1\right)-\left(-h_1\left(t_0\right)\right)-hy\left(t_0\right)d_0^k,-h_1\left(t_2\right)-\left(-h_1\left(t_1\right)\right),\cdots ,-h_1\left(t_N\right)-\left(-h_1\left(t_{N-1}\right)\right)\right)}^{\text{T}}$ ，则上述线性系统(6)和(7)可以写为

$h{B}_N^k{y}_h=H.$ (8)

只要$B_N^k$ 可逆，矩阵$A_{32}$ 就可逆，这只需考虑矩阵$T_N^k$ 和$R_N^k$ 的可逆性。

设$\mathfrak{b}\left(t\right)={\displaystyle \sum _{n=-\infty }^{\infty }{b}_n}{t}^n$ 是Laurent级数，则可定义一个无限维的Toeplitz矩阵

$T\left[\mathfrak{b}\right]=(\begin{array}{cccc}{b}_0& b_{-1}& b_{-2}& \cdots \\ b_1& b_0& b_{-1}& \cdots \\ b_2& b_1& b_0& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}).$

令$T$ 表示复数单位圆，则当t绕$T$ 的逆时针方向移动一圈时，$\mathfrak{b}(t)$ 是一条连续的闭合曲线。令$win{d}_{\mathfrak{b}}$ 表示$\mathfrak{b}\left(t\right)$ 的卷绕数，即曲线$\mathfrak{b}\left(t\right)$ 逆时针绕过原点的次数。特别地，假设对所有的$t\in T$ 有$\mathfrak{b}\left(t\right)\ne 0$ ，且$\mathfrak{b}\left(t\right)$

只有有限多个非零系数，即$\mathfrak{b}(t)={\displaystyle \sum _{j=-r}^s{b}_j{t}^j}$ ，经计算可得$\mathfrak{b}\left(t\right)=t^{-r}{b}_s{\displaystyle \prod _{j=1}^J\left(t-{\delta }_j\right)}{\displaystyle \prod _{i=1}^I\left(t-{\mu }_i\right)}=t^{-r}\overline{\mathfrak{b}}\left(t\right)$ ，其中

$\left|{\delta }_j\right|<1$ ，$\left|{\mu }_i\right|>1$ ，$\overline{\mathfrak{b}}\left(t\right)$ 称为$\mathfrak{b}\left(t\right)$ 的条件多项式。

令${\mathfrak{p}}_j\left(t\right)=t-{\delta }_j$ ，${\mathfrak{q}}_i\left(t\right)=t-{\mu }_i$ ，则$win{d}_{{\mathfrak{p}}_j}=1$ ，$win{d}_{{\mathfrak{q}}_i}=0$ ，$j=1,\cdots ,J$ ，$i=1,\cdots ,I$ ，从而$win{d}_{\mathfrak{b}}=J-r$ 。此外，因为$\mathfrak{b}\left(t\right)\ne 0\left(t\in T\right)$ ，$T\left[\mathfrak{b}\right]$ 是Fredholm算子，即存在算子B，使得$BT\left[\mathfrak{b}\right]-I$ 和$T\left[\mathfrak{b}\right]B-I$ 是紧的[26]，从而可得下面的引理。

引理2.1 [26]算子$T\left[\mathfrak{b}\right]$ 在${\mathfrak{l}}^p\left(1\le p\le \infty \right)$ 上是可逆算子，当且仅当对任意的$t\in T$ ，$\mathfrak{b}(t)\ne 0$ 以及$win{d}_{\mathfrak{b}}=0$ 。

用$T\left[\mathfrak{b}\right]\left(1:n,1:n\right)$ 表示$T\left[\mathfrak{b}\right]$ ，有$\underset{x\to \infty }{\lim }{T}_n\left[\mathfrak{b}\right]=T\left[\mathfrak{b}\right]$ 。因此，$T_n\left[\mathfrak{b}\right]$ 的可逆性由$T\left[\mathfrak{b}\right]$ 来决定。

引理2.2 [26]对Laurent级数$\mathfrak{b}(t)$ ，若$T\left[\mathfrak{b}\right]$ 是可逆算子，则$\underset{n\to \infty }{\lim \sup }\Vert {\left(T_n\left[\mathfrak{b}\right]\right)}^{-1}\Vert <\infty$ ；若$T\left[\mathfrak{b}\right]$ 是不可逆的算子，则$\underset{n\to \infty }{\lim }\Vert {\left(T_n\left[\mathfrak{b}\right]\right)}^{-1}\Vert =\infty$ 。

另一方面，由于基函数线性无关，故$R_N^k$ 是非奇异的。因此，当${\mathfrak{d}}_{{\text{c}}_{N-k-1}^k}\left(t\right)\ne 0$ 且$win{d}_{{\text{c}}_{N-k-1}^k}=0$ 时，线性系统(8)的系数矩阵是可逆的。

引理2.3 [27]假设$f\left(t\right)\in C^d\left[a,b\right],1\le d\le m$ ，对任意的$a\le {\xi }_1<\cdots <{\xi }_m\le b$ ，函数$f\left(t\right)$ 与其关于节点$\left\{{\xi }_i\right\}$ 的$m-1$ 次Lagrange 插值多项式的插值误差记为${\epsilon }_m\left(f;t\right)$ ，即

${\epsilon }_m\left(f;t\right)=f\left(t\right)-{\displaystyle \sum _{j=1}^m{L}_j\left(t\right)f\left({\xi }_j\right),t\in \left[a,b\right]},$

那么${\epsilon }_m\left(f;t\right)$ 可以表示为积分的形式

${\epsilon }_m\left(f;t\right)={\displaystyle {\int }_a^b{\kappa }_d\left(t,s\right)f^d\left(s\right)\text{d}s},t\in \left[a,b\right],$

这里Peano核${\kappa }_d\left(t,s\right)$ 为

${\kappa }_d\left(t,s\right):=\frac{1}{\left(d-1\right)!}\left\{{\left(t-s\right)}_{+}^{d-1}-{\displaystyle \sum _{j=1}^m{L}_j}\left(t\right){\left({\xi }_j-s\right)}_{+}^{d-1}\right\}$ ，其中${\left(t-s\right)}_{+}^p:=\{\begin{array}{l}0,t<s,\hfill \\ {\left(t-s\right)}^p,t\ge s.\hfill \end{array}$

定理2.1 假设在积分代数方程(5)中，核函数都为1，函数$f\left(t\right),g\left(t\right),h_1\left(t\right)\in C^{k+2}\left(\left[0,T\right]\right)$ ，以及${\text{c}}_{N-k-1}^k$ 的卷绕数为0，那么当步长h充分小时，配置边值法存在唯一解，且配置误差为

${\Vert E{1}_N\Vert }_{\infty }=O\left(h^k\right),{\Vert E{2}_N\Vert }_{\infty }=O\left(h^{k+1}\right),{\Vert E{3}_N\Vert }_{\infty }=O\left(h^{k-1}\right),$

其中

$E{1}_N={\left(e_{h1}\left(t_1\right),e_{h1}\left(t_2\right),\cdots ,e_{h1}\left(t_N\right)\right)}^{\text{T}},E{2}_N={\left(e_{h2}\left(t_1\right),e_{h2}\left(t_2\right),\cdots ,e_{h2}\left(t_N\right)\right)}^{\text{T}},$

$E{3}_N={\left(e_{h3}\left(t_1\right),e_{h3}\left(t_2\right),\cdots ,e_{h3}\left(t_N\right)\right)}^{\text{T}}.$

证明：定义误差函数$e_{h1}\left(t\right):=x\left(t\right)-x_h\left(t\right)$ ，$e_{h2}\left(t\right):=y\left(t\right)-y_h\left(t\right)$ ，则有

$0=e_{h2}\left(t\right)+{\displaystyle {\int }_0^t{e}_{h1}\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{e}_{h2}\left(s\right)\text{d}s},$

$0={\displaystyle {\int }_0^t{e}_{h2}\left(s\right)\text{d}s},$

从而可以得到上述误差方程的配置方程

$0=e_{h2}\left(t_n\right)+h{\displaystyle \sum _{j=0}^{n-1}{\displaystyle {\int }_0^1{e}_{h1}\left(t_j+sh\right)\text{d}s}}+h{\displaystyle \sum _{j=0}^{n-1}{\displaystyle {\int }_0^1{e}_{h2}\left(t_j+sh\right)\text{d}s}},$

$0=h{\displaystyle \sum _{j=0}^{n-1}{\displaystyle {\int }_0^1{e}_{h2}\left(t_j+sh\right)\text{d}s}},$

在上式中用n代替n − 1，然后再相减得

$0=e_{h2}\left(t_{n+1}\right)-e_{h2}\left(t_n\right)+h{\displaystyle {\int }_0^1{e}_{h1}\left(t_n+sh\right)\text{d}s}+h{\displaystyle {\int }_0^1{e}_{h2}\left(t_n+sh\right)\text{d}s},$

$0=h{\displaystyle {\int }_0^1{e}_{h2}\left(t_n+sh\right)\text{d}s},$

根据引理2.3，当$n=0,1,\cdots ,N-k-2$ 时，

$0=h\left({\displaystyle {\int }_0^1{\displaystyle \sum _{i=0}^{k+1}{\varphi }_i^k\left(s\right)e_{h2}\left(t_{n+i}\right)}\text{d}s}+O\left(h^{k+2}\right)\right),$

当$n=N-k-1,N-k,\cdots ,N-1$ 时，

$0=h\left({\displaystyle {\int }_0^{n-N+k+2}{\displaystyle \sum _{i=0}^k{\varphi }_i^{k-1}\left(s\right)e_{h2}\left(t_{N-k+i}\right)}\text{d}s}+O\left(h^{k+1}\right)\right),$

因此，对$i=0,1,\cdots ,N-1$ ，有$e_{h2}\left(t_i\right)=O\left(h^{k+1}\right)$ ，进一步可得

$B_N^{k}E{2}_N=\widehat{r},$

其中$\widehat{r}$ 的元素是$O\left(h^{k+1}\right)$ ，由引理2.2知，当$N\to \infty$ 时，${\left(B_N^k\right)}^{-1}$ 的范数是有界的，可得

${\Vert E{2}_N\Vert }_{\infty }=O\left(h^{k+1}\right),h\to 0.$

进一步，我们可以得到，对于$n=0,\cdots ,N-k-2$ ，有$h{\displaystyle {\int }_0^1{\displaystyle \sum _{i=0}^{k+1}{e}_{h1}\left(t_{n+i}\right){\varphi }_i^k\left(s\right)\text{d}s}}=O\left(h^{k+1}\right)$ ，对于$n=N-k-1,N-k,\cdots ,N-1$ ，有$h{\displaystyle {\int }_0^1{\displaystyle \sum _{i=0}^k{e}_{h1}\left(t_{n+i}\right){\varphi }_i^{k-1}\left(s\right)\text{d}s}}=O\left(h^k\right)$ ，从而有

$B_N^{k}E{1}_N={r^{\prime }}_{\widehat{n}},$

其中${r^{\prime }}_{\widehat{n}}$ 的元素是$O\left(h^k\right)$ 。进一步可以得到

${\Vert E{1}_N\Vert }_{\infty }=O\left(h^k\right),h\to 0.$

类似可得

${\Vert E{3}_N\Vert }_{\infty }=O\left(h^{k-1}\right),h\to 0.$

4. 数值实验

下面通过数值例子来验证k−步配置边值方法(k-CBVM)求解方程组(1)的有效性。这里的数值实验都是在Matlab R2018b中实现的，其中积分的计算采用Matlab中自带的积分程序quadgk函数。利用绝对误差的无穷大范数来描述配置边值方法的收敛速度。观察在选取不同步数k时，随着节点数N的增加，数值方法绝对误差的变化情况及其相应的收敛阶。在这里，误差为配置解产生的绝对误差的无穷范数，收敛阶由下式计算：

${\log }_2\frac{erro{r}_{previous}}{erro{r}_{current}}/{\log }_2\frac{N_{current}}{N_{previous}}.$

例1 考虑如下指标-3型积分代数方程

$\{\begin{array}{l}x\left(t\right)=f\left(t\right)+{\displaystyle {\int }_0^t\left(t+s\right)x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t{\left(s+1\right)}^{2}y\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t\left(s^2+t^2+2\right)z\left(s\right)\text{d}s},\\ y\left(t\right)=g\left(t\right)+{\displaystyle {\int }_0^t{\text{e}}^{t+s}x\left(s\right)\text{d}s}+{\displaystyle {\int }_0^t\left(s^2+1\right)y\left(s\right)\text{d}s},\\ 0=h_1\left(t\right)+{\displaystyle {\int }_0^t{\text{e}}^{s+1}y\left(s\right)\text{d}s}\end{array}$

其中$t\in I=\left[0,1\right]$ ，

$f\left(t\right)=t+\left(t^2+2t-1\right)\sin t+2\cos t+\left(2t^2-2t+4\right){\text{e}}^t-t^2-6,$

$g\left(t\right)=-\cos t+\frac{1}{2}\left(\sin t-\cos t\right){\text{e}}^{2t}+\frac{1}{2}{\text{e}}^t+\left(t^2-1\right)\sin t+2t\cos t,$

$h_1\left(t\right)=\frac{1}{2}\left({\text{e}}^{t+1}\cos t+{\text{e}}^{t+1}\sin t-\text{e}\right),$

该方程的解析解为$x\left(t\right)=\sin t$ ，$y\left(t\right)=\cos t$ ，$z\left(t\right)={\text{e}}^t$ 。

Table 1. The variation of the maximum absolute error of x(t) and the convergence order in Example 1

表1. 例1中x(t)的绝对误差的最大值的变化情况以及收敛阶

N	$k=1$		$k=2$		$k=3$
	误差	收敛阶	误差	收敛阶	误差	收敛阶
16	1.50 × 10^−3	—	1.47 × 10^−4	—	1.10 × 10^−5	—
32	6.20 × 10^−4	1.27	3.98 × 10^−5	1.89	1.19 × 10^−6	3.21
48	3.87 × 10^−4	1.17	1.81 × 10^−5	1.94	3.34 × 10^−7	3.14
64	2.80 × 10^−4	1.12	1.03 × 10^−5	1.96	1.37 × 10^−7	3.10
80	2.19 × 10^−4	1.10	6.64 × 10^−6	1.97	6.88 × 10^−8	3.08
96	1.80 × 10^−4	1.08	4.63 × 10^−6	1.98	3.93 × 10^−8	3.07

Table 2. The variation of the maximum absolute error of y(t) and the convergence order in Example 1

表2. 例1中y(t)的绝对误差的最大值的变化情况以及收敛阶

N	$k=1$		$k=2$		$k=3$
	误差	收敛阶	误差	收敛阶	误差	收敛阶
16	3.21 × 10^−4	—	2.95 × 10^−5	—	1.26 × 10^−6	—
32	7.95 × 10^−5	2.01	3.78 × 10^−6	2.96	7.48 × 10^−8	4.07
48	3.52 × 10^−5	2.01	1.13 × 10^−6	2.98	1.45 × 10^−8	4.04
64	1.98 × 10^−5	2.01	4.79 × 10^−7	2.99	4.55 × 10^−9	4.03
80	1.27 × 10^−5	2.00	2.46 × 10^−7	2.99	1.85 × 10^−9	4.02
96	8.78 × 10^−6	2.00	1.42 × 10^−7	2.99	8.90 × 10^−10	4.02

Table 3. The variation of the maximum absolute error of z(t) and the convergence order in Example 1

表3. 例1中z(t)的绝对误差的最大值的变化情况以及收敛阶

N	$k=1$		$k=2$		$k=3$
	误差	收敛阶	误差	收敛阶	误差	收敛阶
16	5.50 × 10^−^3	—	8.08 × 10^−4	—	9.66 × 10^−5	—
32	5.60 × 10^−3	−0.03	4.52 × 10^−4	0.84	1.97 × 10^−5	2.29
48	5.40 × 10^−3	0.09	3.15 × 10^−4	0.89	8.14 × 10^−6	2.18
64	5.30 × 10^−3	0.07	2.42 × 10^−4	0.92	4.41 × 10^−6	2.13
80	5.30 × 10^−3	0.00	1965 × 10^−4	0.94	2.76 × 10^−6	2.10
96	5.20 × 10^−3	0.10	1.65 × 10^−4	0.95	1.88 × 10^−6	2.09

表1、表2和表3分别给出了例1中$x\left(t\right),y\left(t\right)$ 和$z\left(t\right)$ 在取不同步数k下，绝对误差的最大值随节点数N的变化情况以及相应收敛阶的变化情况。通过观察表1、表2、表3中的数据，可以发现，当步数k不变时，随着节点数N的增加，相应的误差在逐渐减小；当节点数N不变时，随着步数k的增加，相应的误差也在在逐渐减小。并且随着节点数N的增加，三个表中的收敛阶分别靠近$k,k+1$ 和$k-1$ 。由此可知，利用配置边值方法求解指标-3型积分代数方程是有效的，并且具有较快的收敛速度。

5. 结论

本文通过对区间采用网格划分，基于多步配置方法，提出了一类数值求解指标-3型积分代数方程的k−步配置边值方法。首先利用Lagrange基函数构造配置边值法，将原方程离散成一个线性系统，后通过整体求解线性系统得到配置解。然后研究了该方法求解指标-3型积分代数方程的可解性和收敛性。最后利用Matlab进行数值实验，实验结果表明，运用此类k−步配置边值方法求解指标-3型积分代数方程是有效的，并且可以达到较高的收敛阶。

基金项目

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