Gauss型积分公式

摘要求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。这时我们可以通过数值方法求出函数积分的近似值。当然再用近似值代替真实值时，误差精度是我们需要考虑因素，但是除了误差精度以外，还可以用代数精度来判断其精度的高低。已知n+1点的Newton-Cotes型积分公式，当n为奇数时，其代数精度为n；当n为偶数时，其代数精度达到n+1。若对随机选取的n+1个节点作插值型积分公式也仅有n次代数精度。如何选取适当的节点，能使代数精度提高？Gauss型积分公式可是实现这一点，但是Gauss型求积公式，需要被积函数满足的条件是正交，这一条件比较苛刻。因此本实验将针对三种常用的Gauss型积分公式进行讨论并编程实现。关键词：Newton-Cotes型积分公式正交多项式代数精度11、实验目的1)通过本次实验体会并学习Gauss型积分公式，在解决如何取节点能提高代数精度这一问题中的思想方法。2)通过对Gauss型积分公式的三种常见类型进行编程实现，提高自己的编程能力。3)用实验报告的形式展现，提高自己在写论文方面的能力。2、算法流程下面介绍三种常见的Gauss型积分公式1)高斯-勒让德（Gauss-Legendre）积分公式勒让德（Legendre）多项式如下定义的多项式Ln(x)=12n𝑛!𝑑𝑛𝑑𝑥𝑛(𝑥2−1)𝑛,𝑥∈[−1,1],𝑛=0,1,2⋯称作勒让德多项式。由于(𝑥2−1)𝑛是2n次多项式，所以Ln(x)是n次多项式，其最高次幂的系数An与多项式12n𝑛!𝑑𝑛𝑑𝑥𝑛(𝑥(2n))=12n𝑛!2𝑛(2𝑛−1)(2𝑛−2)⋯(𝑛+1)𝑥𝑛的系数相同。也就是说n次勒让德多项式具有正交性即勒让德多项式Ln(x)是在[−1,1]上带𝜌(x)=1的n次正交多项式，而且(𝐿𝑚,𝐿𝑛)=∫𝐿𝑚(𝑥)𝐿𝑛(𝑥)𝑑𝑥1−1={0,𝑚≠𝑛22𝑛+1,𝑚=𝑛这时Gauss型积分公式的节点就取为上述多项式Ln(x)的零点，相应的Gauss型积分公式为∫𝑓(𝑥)𝑑𝑥1−1≈∑𝐴𝑘𝑓(𝑥𝑘)𝑛𝑘=1此积分公式即成为高斯-勒让德积分公式。其中Gauss-Legendre求积公式的系数𝐴𝑘=∫𝜌(𝑥)𝜔𝑛(𝑥)(𝑥−𝑥𝑘)𝜔′𝑛(𝑥)𝑑𝑥1−1=∫𝜌(𝑥)𝐿𝑛(𝑥)(𝑥−𝑥𝑘)𝐿′𝑛(𝑥)𝑑𝑥1−12其中k的取值范围为𝑘=1,2,⋯,𝑛Gauss点和系数不容易计算，但是在实际计算中精度要求不是很高，所以给出如下表所示的部分Gauss点{𝑥𝑘}和系数{𝐴𝑘}，在实际应用中只需查表即可。nxAnxA1026±0.9324695142±0.6612093865±1.23861918160.1713244920.3607615730.4679139342±0.577350269217±0.9491079123±0.7415311856±0.405845151400.1294849660.2797053910.3818300500.4179591833±0.7745966692000.55555555560.88888888894±0.8611363116±0.33998104360.34785484510.65214515498±0.9602898565±0.7966664774±0.5255324099±0.18343464250.1012285360.2223810340.3137066450.3626837835±0.9061798459±0.538469310100.23692688510.47862867050.56888888892)高斯-拉盖尔（Gauss-Laguerre）积分公式拉盖尔（Laguere）多项式Ln(x)=e𝑥𝑑𝑛𝑑𝑥𝑛(𝑥𝑛𝑒−𝑥),0≤𝑥+∞,𝑛=0,1,2⋯称为拉盖尔多项式。其首项系数为(−1)𝑛,且具有性质：正交性，在区间[0,+∞)上关于权函数𝜌(x)=e−𝑥正交，而且(𝐿𝑚,𝐿𝑛)=∫e−𝑥𝐿𝑚(𝑥)𝐿𝑛(𝑥)𝑑𝑥∞0={0,𝑚≠𝑛(𝑛!)2,𝑚=𝑛积分区间为[0,+∞)，权函数为𝜌(x)=e−𝑥的Gauss型积分公式称为高斯-拉盖尔积分公式，其中Gauss点为拉盖尔多项式Ln(x)的零点，高斯-拉盖尔积分公式为∫e−𝑥𝑓(𝑥)𝑑𝑥∞0≈∑𝐴𝑘𝑓(𝑥𝑘)𝑛𝑘=1同样高斯-拉盖尔积分公式的Gauss点和求积系数如下表所示：3nxAnxA20.58578643763.41421356240.85355339050.146446609450.26356031971.41340305913.59642577107.085810005812.64080084420.52175561050.39866681100.07594244970.00361175870.000023370030.41577455672.29428036026.28994508290.71109300990.278517733540.32254768961.74576110114.53662029699.39507091230.60315410430.35741869240.03888790850.000539294760.22284660411.18893210162.99273632605.77514356919.837467418315.98287398060.45896467930.41700083070.11337338200.01039919750.00026101720.00000089853)高斯-埃尔米特（Gauss-Hermite）积分公式埃尔米特（Hermite）多项式Hn(x)=(−1)𝑛e𝑥2𝑑𝑛e−𝑥2𝑑𝑥𝑛,−∞𝑥+∞,𝑛=0,1,2⋯被称作埃尔米特多项式，其首项系数为2𝑛，具有性质如下正交性，在区间(−∞,+∞)上关于权函数e−𝑥2正交，而且(𝐻𝑚,𝐻𝑛)=∫e−𝑥2𝐻𝑚(𝑥)𝐻𝑛(𝑥)𝑑𝑥+∞−∞={0,𝑚≠𝑛2𝑛𝑛!√𝜋,𝑚=𝑛积分区间为(−∞,+∞)，权函数为𝜌(x)=e−𝑥2的Gauss型积分公式称为Gauss-Hermite积分公式，其Gauss点就是Hermite正交多项式Hn(x)的零点。Gauss-Hermite求积公式为∫e−𝑥2𝑓(𝑥)𝑑𝑥∞−∞≈∑𝐴𝑘𝑓(𝑥𝑘)𝑛𝑘=1同样高斯-埃尔米特积分公式的Gauss点和求积系数如下表所示：4nxAnxA2±0.70710678110.88622692557±0.4360774119±1.3358490704±2.35030497360.72462959520.15706732030.00453000993±1.224744871400.29540897511.81635900064±0.5246476232±1.65068012380.80491409000.08131283548±0.8162878828±1.6735516287±2.651961356300.42560725260.05451558280.00097178120.81026461755±0.9585724646±2.02018287040.39361932310.01995324210.94530872040.56888888893、算法实例1)用3点Gauss型求积公式计算∫cosx𝑑𝑥1−1解：根据积分限可以知道应该用Gauss-Legendre积分公式，具体程序如下所示#includeiostream#includemath.husingnamespacestd;constintM(10);voidmain(){inti=0;intn=0;intm=0;intsign=0;doublesum=0;doublex[M]={0};doubleA[M]={0};doublex1[]={0};doublex2[]={-0.57735502692,0.57735502692};doublex3[]={-0.77459666920,0.77459666920,0};doublex4[]={-0.8611363116,0.8611363116,-0.3399810436,0.33998104365};doublex5[]={-0.9061798459,0.9061798459,-0.53846931010,0.53846931010,0};doublex6[]={-0.9324695142,0.9324695142,-0.6612093865,0.6612093865,-1.2386191816,1.2386191816};doublex7[]={-0.9491079123,0.9491079123,-0.7415311856,0.7415311856,-0.40584515140,0.40584515140,0};doublex8[]={-0.9602898565,0.9602898565,-0.7966664774,0.7966664774,-0.5255324099,0.5255324099,-0.1834346425,0.1834346425};doubleA1[]={2};doubleA2[]={1};doubleA3[]={0.5555555556,0.8888888889};doubleA4[]={0.3478548451,0.6521451549};doubleA5[]={0.2369268851,0.4786286705,0.5688888889};doubleA6[]={0.1713244924,0.3607615730,0.4679139346};doubleA7[]={0.1294849662,0.2797053915,0.3818300505,0.4179591834};doubleA8[]={0.1012285363,0.2223810345,0.3137066459,0.3626837834};cout请输入节点个数endl;cinn;switch(n){case1:for(i=0;in;i++){x[i]=x1[i];A[i]=A1[i];}break;case2:for(i=0;in;i++){x[i]=x2[i];}for(i=0;in;i++){A[i]=A2[i/2];}break;case3:for(i=0;in;i++){x[i]=x3[i];}for(i=0;in;i++){A[i]=A3[i/2];}break;6case4:for(i=0;in;i++){x[i]=x4[i];}for(i=0;in;i++){A[i]=A4[i/2];}break;case5:for(i=0;in;i++){x[i]=x5[i];}for(i=0;in;i++){A[i]=A5[i/2];}break;case6:for(i=0;in;i++){x[i]=x6[i];}for(i=0;in;i++){A[i]=A6[i/2];}break;case7:for(i=0;in;i++){x[i]=x7[i];}for(i=0;in;i++){A[i]=A7[i/2];}break;case8:for(i=0;in;i++){x[i]=x8[i];}for(i=0;in;i++){A[i]=A8[i/2];}break;default:cout输入出错，请从新输入！！endl;break;}for(i=0