Cauchy核奇异积分的反Gauss求积算法

李寒嫣，张彦铎

武汉工程大学计算机科学与工程学院，湖北 武汉 430205

近年来，奇异积分理论和它的数值计算是国际数学逼近论界比较活跃的研究课题。奇异积分已经在许多领域得到了广泛的应用，包括信号处理、计算机辅助工程（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求积公式

理论已广泛应用于工程技术计算和数值理论研究中，其中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］

其中xk仍是式（1）中求积结点，τj是Stieltjes正交多项式的零点，求积公式（2）的代数精度为3n+1，用近似估计式（1）中的求积误差。在实际计算中，给误差估计带来很多方便。Gauss-Kronrod求积公式（2）也有不足［11-12］，求积权σk不一定都为正，对某些权函数求积公式不一定存在。1996年，Laurie首次提出正常积分内插型反Gauss求积公式［13］

结点tk，k=1,2,…,n+1是一正交多项式的零点，且-1＜t1＜x1＜t2＜x2＜…＜t n＜xn＜t n+1＜1，ωk＞0。因不涉及另外的正交多项式，反Gauss求积公式比Gauss-Kronrod求积公式更容易构造，式（3）中误差与式（1）中误差反号，即

1.2 正交多项式

首一Legendre正交多项式满足三项递推关系

Laurie定义新的多项式［13］，

定理1［13-15］正交多项式πAGn+1(t)的零点tk，k=1,2,…,n+1都为实结点且与Legendre多项式πn(x)的零点xi，i=1,2,…,n交错分布，且

2 Cauchy核奇异积分的反Gauss公式求积

为保证含Cauchy核奇异积分

存在，要求密度函数f满足Hölder条件，记作f∈H。设f∈Ck+1[-1,1]，令

于是Fλ∈Ck[-1,1]，当f为一n次多项式时，Fλ(t)是n-1次多项式。利用分离奇异性，I(f,λ)=-1＜λ＜1。当f∈H时是一正常广义积分，利用正常积分反Gauss求积公式可构造I(Fλ)的反Gauss求积公式，从而得到I(f,λ)的反Gauss求积公式。 由式（3）其 中当λ=ti0时 ，

于是有下面结果。

定理2设f∈C2n+1[-1,1]，I(f,λ)的反Gauss求积公式：

其中已取定logz的某一连续分支，结点t1＜t2＜…＜tn+1为多项式的零点，

证设函数1,2,...,n+1，其中tμ为正交多项式（6）的零点，由奇异积分定义，

由式（4）得

定理3当定理2中条件满足时，I(f,λ)的余项积分表达式为

证设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插值多项式，于是

式中右侧第一个积分对2n-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软件进行仿真实验，通过数值实例的误差数据表和逼近图像，说明求积公式的误差渐近性。误差定义为其中I(f,λ)表示积分精确值。实验流程如图1所示。

图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函数，奇异点t=0.5，考虑积分

取求积结点个数N为31,35,…,63，实验结果如图2所示。

图2 积分渐近图Fig.2 Integral asymptotic graph

由图2知，求积公式计算的积分值随着求积结点数目增多而不断逼近积分的精确值，误差随求积结点数增大而逐渐减小；计算误差与求积结点的个数成反比。同时观察与奇异积分Gauss求积算法所得的误差符号相反，数值相等。这些与理论上求积公式的误差渐近分析一致。

图3为奇异积分的Gauss求积与反Gauss求积的n+1个点的误差。误差曲线接近，但本文所论述的奇异积分反Gauss算法在计算过程上更为方便，说明建立的求积公式具有可行性。

图3 n+1个结点的积分误差Fig.3 Integral errors of n+1 nodes

4 结 论

本文创建的奇异积分反Gauss求积算法在工程技术计算中有良好应用前景，为计算机某些工程应用软件的开发提供了理论依据。对奇异积分进行数值求积时，随着求积节点数目N的增多，通过反Gauss求积公式（10）运算后得到的误差也在缩小，误差曲线也较为平滑，所得积分近似值逐渐逼近积分精确值。但是由于实验设备运行能力不足，求积公式的收敛速度本文暂未涉及，且选取的求积节点数较少，所得的一些数据略显粗糙，如果选取更多的节点，得到的实验结果将会更好。