近年来，奇异积分理论和它的数值计算是国际数学逼近论界比较活跃的研究课题。奇异积分已经在许多领域得到了广泛的应用，包括信号处理、计算机辅助工程（computer aided engineering，CAE）等［1-3］。在奇异积分理论被广泛研究的同时，它的数值计算也随之被提了出来［4］。国内比较有代表性的是武汉大学路见可教授和杜金元教授在这方面的研究工作［5-6］。近十年来，由于反Gauss求积算法具有良好性质，其应用范围更广，对它的研究十分活跃［7-8］。目前对奇异积分的反Gauss求积算法研究文献尚少， 2018年文献［9］用分离奇异性的方法，讨论了含Cauchy核奇异积分的Gauss-Kronrod求积问题。本文讨论含Cauchy核奇异积分的反Gauss求积算法，建立奇异积分的反Gauss求积公式，讨论了求积公式的代数精度和余项积分表达式等问题，在计算机上实现数值仿真。
1 相关研究
1.1 Gauss型求积公式
正常积分Gauss求积公式
[I(f)=-1 1f(t)dt=k=1nAkf(xk)+RGn(f)=] [QGn(f)+RGn(f)] （1）
理论已广泛应用于工程技术计算和数值理论研究中，其中[x1, x2, ?, xn]为[n]次Legendre多项式的零点，求积权[Ak>0]，式（1）的代数精度为[2n-1]。为了估计误差[RGn(f)]，通常采用[RGn(f)≈|QGm(f)-]
[QGn(f)|]，[m>>n]，这种近似替代需要重新计算[f]在[m]个结点的值及求积权，这将增加更大的计算量。1964年，Kronrod为了避免重复计算[f(x)]在[m]个结点上的值及求积权，首次构造出Gauss-Kronrod求积公式［10］
[-1 1f(t)dt=k=1nAkf(xk)+j=1n+1σjf(τj)+]
[RKn(f)=QKn(f)+RKn(f)] （2）
其中[xk]仍是式（1）中求积结点，[τj]是Stieltjes正交多项式的零点，求积公式（2）的代数精度为[3n+1]，用[RGn(f)≈QGn(f)-QKn(f)]近似估计式（1）中的求积误差。在实际计算中，给误差估计带来很多方便。Gauss-Kronrod求积公式（2）也有不足［11-12］，求积权[σk]不一定都为正，对某些权函数求积公式不一定存在。1996年，Laurie首次提出正常积分内插型反Gauss求积公式［13］
[-1 1f(t)dt=k=1n+1ωkf(tk)+RAGn+1(f)=]
[QAGn+1(f)+RAGn+1(f)] （3）
结点[tk]，[k=1, 2, ?, n+1]是一正交多项式的零点，且-1[<t1<x1<t2<x2<?<tn<xn<tn+1<1]，[ωk>0]。因不涉及另外的正交多项式，反Gauss求积公式比Gauss-Kronrod求积公式更容易构造，式（3）中误差与式（1）中误差反号，即
[I(p)-QAGn+1(p)=-(I(p)-QGn(p))]， [p∈P2n+1]
（4）
1.2 正交多项式
首一Legendre正交多项式满足三项递推关系
[πj+1(t)=tπj(t)-βjπj-1(t)]，[j=0,1,2,…n]，
[π0(t)=1]，[π-1(t)=0]， （5）
其中[βj=I(π2j)I(π2j-1)]，[j=0, 1, ?, n]。
Laurie定义新的多项式［13］，
[πAGn+1(t)=πn+1(t)-βnπn-1(t)] （6）
证明[πAGn+1(t)]是满足正交关系
[-1 1πAGn+1(t)p(t)dt=]
[-1 1[πn+1(t)-βnπn-1(t)]p(t)dt=0]，[p∈Pn-2] （7）
的正交多项式，且有如下定理：
定理1［13-15］ 正交多项式[πAGn+1(t)]的零点[tk]，[k=1,2,…,n+1]都为实结点且与Legendre多项式[πn(x)]的零点[xi]，[i=1,2,…,n]交错分布，且
[-1<t1<x1<t2<x2<?<xn<tn+1<1] （8）
2 Cauchy核奇异积分的反Gauss公式求积
为保证含Cauchy核奇异积分
[I(f, λ)=-1 1f(t)t-λdt=limε→0+-1λ-εf(t)t-λdt+λ+ε 1f(t)t-λdt,]
[-1<λ<1]
存在，要求密度函数[f]满足H?lder条件，记作[f∈H]。设[f∈Ck+1[-1,1]]，令
[Fλ(t)=f(t)-f(λ)t-λ,t≠λ,f(λ),t=λ,]
于是[Fλ∈Ck[-1,1]]，当[f]为一[n]次多项式时，[Fλ(t)]是[n-1]次多项式。利用分离奇异性，[I(f, λ)=]
[-1 1f(t)t-λdt=][-1 1f(t)-f(λ)t-λdt+f(λ)][-1 11t-λdt]，[-1<λ<1]。当[f∈H]时，[I(Fλ)=-1 1f(t)-f(λ)t-λdt]是一正常广义积分，利用正常积分反Gauss求积公式可构造[I(Fλ)]的反Gauss求积公式，从而得到[I(f, λ)]的反Gauss求积公式。 由式（3）[-1 1f(t)-f(λ)t-λdt≈][QAGn+1(f)]，其中[QAGn+1(f)=]
[i=1n+1ωif(ti)-f(λ)ti-λ]，[λ≠ti]，当[λ=ti0]时，[QAGn+1(f)=][i=1i≠i0n+1ωif(ti)-f(λ)ti-λ+][ωi0f(ti0)]。
于是有下面结果。
定理2 设[f∈C2n+1[-1, 1]]，[I(f, λ)]的反Gauss求积公式：
[I(f, λ)=I(Fλ)+f(λ)-1 11t-λdt= ][QAGSn+1(f, λ)+RAGSn+1(f, λ)]，[-1<λ<1] （9）
其中
[QAGSn+1(f, λ)=]
[i=1n+1ωif(ti)-f(λ)ti-λ+f(λ)log1-λ1+λ,λ≠tii=1i≠i0n+1ωif(ti)-f(λ)ti-λ+ωi0f(λ)+f(λ)log1-λ1+λ,λ=ti0 ]
（10）
已取定[logz]的某一连续分支，结点[t1<t2<?<tn+1]为多项式[πAGn+1(t)]的零点，
[ωi=2-1 1π2n(t)dtπn(ti)(πAGn+1)(ti)]， [i=1, 2, ?, n+1]
[QAGSn+1(f, λ)]的代数精度至少为[2n]。
证 设函数[p(t)=πAGn+1(t)t-tμπn(t)]，[μ=1, 2, ..., n+1]，其中[tμ]为正交多项式（6）的零点，由奇异积分定义，
[-1 1πAGn+1(t)t-tμπn(t)dt=]
[ωμ(πAGn+1)′(tμ)πn(tμ)+RAGn+1πAGn+1(t)t-tμπn(t)] （11）
由于[πAGn+1(t)t-tμ]是首一[n]次多项式，[πAGn+1(t)t-tμ=πn(t)+cn-1πn-1(t)+?+c1π1(t)+c0]，其中[c0, c1,?, cn-1]均为常数。 由正交性
[-1 1πAGn+1(t)t-tμπn(t)dt=-1 1π2n(t)dt] （12）
由式（4）得
[RAGn+1πAGn+1(t)t-tμπn(t)][=-RnπAGn+1(t)t-tμπn(t)][=--1 1πAGn+1(t)t-tμπn(t)dt=--1 1π2n(t)dt]， （13）
将式（12）、式（13） 代入式（11），由[-11π2n(t)dt=ωμ(πAGn+1)′(ti)πn(ti)--11π2n(t)dt]得求积系数[ωi]。 [?f∈P2n]，有[f(t)-f(λ)t-λ∈P2n-1]。 根据式（13）得到[QAGSn+1(f, λ)]的代数精度至少为[2n]。
定理3 当定理2中条件满足时，[I(f,λ)]的余项积分表达式为
[RAGSn+1=1(2n)!-1 1F2nλ(ξ)(πAGn+1)2(t)dt]， [-1<ξ<1]，
证 设[H2n-1(t)]为[Fλ(t)]在结点[ti]，[i=1,2,?n+1]满足条件[H2n-1(ti)=Fλ(ti)]，[H′2n-1(ti)=F′λ(ti)]的[2n-1]次Hermite插值多项式，于是
[Fλ(t)=H2n-1(t)+F(2n)λ(ξ)(2n)!(πAGn+1)2(t)]， [-1<ξ<1]，
[-1 1Fλ(t)dt=-1 1H2n-1(t)dt+-1 1F(2n)λ(ξ)(2n)!(πAGn+1)2(t)dt]
式中右侧第一个积分对[2n-1]次多项式精准成立，而
[-1 1F2nλ(ξ)(2n)!(πAGn+1)2(t)dt=][-1 1Fλ(t)dt--1 1H2n-1(t)dt=][-1 1Fλ(t)dt-i=1n+1ωiFλ(ti)]，
于是
[RAGSn+1(f, λ)=I(f,λ)-QAGSn+1(f, λ)=]
[-1 1Fλ(t)dt-i=1n+1ωiFλ(ti)=]
[1(2n)!-1 1F2nλ(ξ)(πAGn+1)2(t)dt]，[-1<ξ<1]。
3 实验部分
3.1 实验配置
联想IdeaPad笔记本电脑一台，该笔记本电脑包含一台英特尔i5-6300HQ物理CPU，该CPU包含4个内核，主频2.30 GHz，内存8 GB，硬盘1 TB，物理网卡1个。安装win10专业版操作系统，使用MATLAB软件对含Cauchy核奇异积分的反Gauss求积公式进行实例数值实现。
3.2 数值仿真实验
对含Cauchy核奇异积分的反Gauss求积公式（9），用MATLAB软件进行仿真实验，通过数值实例的误差数据表和逼近图像，说明求积公式的误差渐近性。误差定义为：[Er=|I(f, λ)-QAGSn+1][f, λ|]，其中[I(f, λ)]表示积分精确值。实验流程如图1所示。
<G:\武汉工程大学\2022\第2期\李寒嫣-1.tif>
图1 实验流程图
Fig. 1 Experimental flowchart
利用MATLAB的自适应性数值求积算法得到公式（9）的仿真，取[ε=10-8]。用t表示奇异点，f表示所求函数，在处理奇异点时的伪代码如下：
sym x；
g=f/（x-t）；
exact=int（subs（g，x），x，-1，t-epsilon）+int（subs（g，x），x，t+epsilon）；
例1 函数[f(x)=x+1]，奇异点[t=0.5]，考虑积分[-11x+1x-0.5dx]。
取求积结点个数[N]为[31,35,?,63]，实验结果如图2所示。
由图2知，求积公式计算的积分值随着求积结点数目增多而不断逼近积分的精确值，误差随求积结点数增大而逐渐减小；计算误差与求积结点的个数成反比。同时观察与奇异积分Gauss求积算法所得的误差符号相反，数值相等。这些与理论上求积公式的误差渐近分析一致。
图3为奇异积分的Gauss求积与反Gauss求积的[n+1]个点的误差。误差曲线接近，但本文所论述的奇异积分反Gauss算法在计算过程上更为方便，说明建立的求积公式具有可行性。
<G:\武汉工程大学\2022\第2期\李寒嫣-3.tif>[0.000 05
0.000 04
0.000 03
0.000 02
0.000 01
0][误差值][反Gauss
Gauss][30 35 40 45 50 55 60 65
节点数]
图3 [n+1]个结点的积分误差
Fig. 3 Integral errors of [n+1] nodes
4 结 论
本文创建的奇异积分反Gauss求积算法在工程技术计算中有良好应用前景，为计算机某些工程应用软件的开发提供了理论依据。对奇异积分进行数值求积时，随着求积节点数目[N]的增多，通过反Gauss求积公式（10）运算后得到的误差也在缩小，误差曲线也较为平滑，所得积分近似值逐渐逼近积分精确值。但是由于实验设备运行能力不足，求积公式的收敛速度本文暂未涉及，且选取的求积节点数较少，所得的一些数据略显粗糙，如果选取更多的节点，得到的实验结果将会更好。