高斯—勒让德积分公式

　摘要：

　高斯—勒让德积分公式可以用较少节点数得到高精度的计算结果,是现在现实生活中经常运用到的数值积分法。然而，当积分区间较大时，积分精度并不理想。

　 The adva ntage of Gauss-Legendre integral formula is tend to get high-precision calculational result by using fewer Gauss-points, real life is now often applied numerical integration method. But the precision is not good when the length of integral interval is longer.

　关键字：

　积分计算，积分公式，高斯—勒让德积分公式，MATLAB

　Keyword:

　Integral Calculation , Integral formula ,Gauss-Legendre integral formula, Matlab

　 引言：

　 众所周知，微积分的两大部分是微分与积分。微分实际上是求一函数的导数，而积分是已知一函数的导数，求这一函数。所以，微分与积分互为逆运算。

　实际上，积分还可以分为两部分。第一种，是单纯的积分，也就是已知导数求原函数，称为不定积分。

　相对而言，另一种就是定积分了，之所以称其为定积分，是因为它积分后得出的值是确定的，是一个数，而不是一个函数。

　计算定积分的方法很多，而高斯—勒让德公式就是其中之一。

　高斯积分法是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到。

　高斯—勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)以此取i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n-1，而且是最高的。通常运用的是（-1，1）的积分节点和积分系数，其他积分域是通过变换x=(b-a)t/2   +(a+b)/2 变换到-1到1之间积分。

　 1. 现有的方法和理论

　1.1高斯 勒让德求积公式 在高斯求积公式(4.5.1)中，若取权函数，区间为，则得公式

　 我们知道勒让德多项式是区间上的正交多项式，因此，勒让德多项式的零点就是求积公式(上式)的高斯点．形如(上式)的高斯公式特别地称为高斯－勒让德求积公式． 若取的零点做节点构造求积公式

　 令它对准确成立，即可定出．这样构造出的一点高斯－勒让德求积公式是中矩形公式．再取的两个零点构造求积公式

　 令它对都准确成立，有 ． 由此解出，从而得到两点高斯－勒让德求积公式 ． 三点高斯－勒让德求积公式的形式是 ． 如表列出高斯－勒让德求积公式的节点和系数．

　0 0.0000000 2.0000000 1 0.5773503 1.0000000 2 0.7745967 0.0000000 0.5555556 0.8888889 3 0.8611363 0.3399810 0.3478548 0.6521452 4 0.9061798 0.5384693 0.0000000 0.2369269 0.4786287 0.5688889 公式(4.5.9)的余项由(4.5.8)得 ， 这里是最高项系数为１的勒让德多项式，由(3.2.6)及(3.2.7)得 ．　　　

　当时，有 ． 它比辛普森公式余项还小，且比辛普森公式少算一个函数值． 当积分区间不是[－１，１]，而是一般的区间时，只要做变换

　可将化为［－１，１］，这时 ．　

　 对等式右端的积分即可使用高斯－勒让德求积公式．

　1.2复化Gauss-Legendre求积公式 将被积区间m等分, 记, 作变换

　在每个小区间上应用Gauss-Legendre公式, 累加即得复化Gauss-Legendre求积公式

　不妨设

　则有: Gauss点个数时,

　Gauss点个数时,

　总结复化Gauss-Legendre求积过程如下: 1. 分割区间, 记录区间端点值； 2. 通过查表或求解非线性方程组, 在所有小区间上, 将Gauss系数和Gauss点的值代入变量替换后的公式； 3. 将所有区间的结果累加, 即得到整个区间上的积分近似值.

　针对Gauss点个数和的复化Gauss-Legendre求积公式编写的一个简单的MATLAB函数 compgauss() 如下:

　function [ ] = compgauss(a, b, n) % Composite Gauss Integration % Equation Type: n=2, n=3 % Coded by Nan.Xiao 2010-05-25 % Step.1 Divide Interval % Step.2 Calculate % Step.3 Sum Results format long f = @(x) exp(x).*sin(x); h=(b-a)/n; xk=zeros(n+1,1); xk(1,1)=a; xk(n+1,1)=b; fk1=zeros(n,1); fk2=zeros(n,1); for i=1:n-1

　 xk(i+1,1)=a+h*i; end for j=1:n

　 fk1(j)=f((xk(j)+xk(j+1))/2+(h/2)*(-1/sqrt(3)))+...

　 f((xk(j)+xk(j+1))/2+(h/2)*(1/sqrt(3))); end for r=1:n

　 fk2(r)=(5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(-sqrt(15)/5))+...

　 (8/9)*f((xk(r)+xk(r+1))/2+(h/2)*(0))+...

　 (5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(sqrt(15)/5)); end mysum1=h*sum(fk1)/2; mysum2=h*sum(fk2)/2; disp(＇Result of 2 Nodes:＇) disp(mysum1); disp(＇Result of 3 Nodes:＇) disp(mysum2); end

　1.3龙贝格，三点，五点以及变步长高斯勒让德求积法 以下是关于龙贝格，三点，五点以及变步长高斯勒让德之间精度的相互比较 #include <iostream.h>

　#include <math.h>

　#include <iomanip.h>

　#define Precision1 0.000000000001

　# define e         2.71828183

　#define  MAXRepeat 10  

　double function (double x)

　{

　  double s;

　s=1/x;

　return s;

　}

　double Romberg(double a,double b,double f(double x))

　{

　  int m,n,k;

　 double y[MAXRepeat],h,ep,p,xk,s,q;

　h=b-a;

　 y[0]=h*(f(a)+f(b))/2.0;//计算T`1`(h)=1/2(b-a)(f(a)+f(b));

　 m=1;

　n=1;

　ep=Precision1+1;

　 while((ep>=Precision1)&&(m<MAXRepeat))

　{

　 p=0.0;

　for(k=0;k<n;k++)

　{

　 xk=a+(k+0.5)*h; 

　p=p+f(xk);  

　        }       

　 p=(y[0]+h*p)/2.0;  //T`m`(h/2),变步长梯形求积公式

　 s=1.0;

　for(k=1;k<=m;k++)

　{

　s=4.0*s;// pow(4,m)

　q=(s*p-y[k-1])/(s-1.0);

　y[k-1]=p;

　p=q;

　}

　 ep=fabs(q-y[m-1]);

　 m=m+1;            

　y[m-1]=q;

　n=n+n;   //  2 4 8 16 

　  h=h/2.0;//二倍分割区间

　return q;

　}

　double ThreePointGaussLegendre(double a,double b,double f(double x))

　{

　 double x,w;

　 static double X[3]={-sqrt(15)/5.0,0,sqrt(15)/5.0};

　  static double L[3]={5/9.0,8/9.0,5/9.0};

　    w=0.0;

　  for(int i=0;i<3;i++)

　          {

　             x=((b-a)*X[i]+(b+a))/2.0;               

　           w=w+f(x)*L[i];

　          }

　    return w;

　}

　 double FivePointGaussLegendre(double a,double b,double f(double x))

　{

　  double x,w;

　    static double X[5]={-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459};

　    static double L[5]={0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851};

　  w=0.0;

　    for(int i=0;i<5;i++)

　          {

　             x=((b-a)*X[i]+(b+a))/2.0;               

　              w=w+f(x)*L[i];//每一次小区间利用勒让德公式计算的结果

　        }

　   return w;

　}

　double FivePointPrecisionGaussLegendre(double a,double b,double f(double x))

　{

　  int m,i,j;

　    double s,p,ep,h,aa,bb,w,x,g;

　    static double X[5]={-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459};

　   m=1;

　  h=b-a;

　    s=fabs(0.001*h);

　p=1.0e+35;     ep=Precision1+1;

　   while((ep>=Precision1)&&(fabs(h)>s))     {

　g=0.0;

　        for(i=0;i<m;i++)

　        {

　           aa=a+i*h;             bb=aa+h;

　            w=0.0;

　         for(j=0;j<=4;j++)

　           {

　              x=((bb-aa)*X[j]+(bb+aa))/2.0;                 w=w+f(x)*L[j];

　           }

　          g=g+w;//各个区间计算结果之和相加

　        }

　           g=g*h/2.0;

　            ep=fabs(g-p)/(1.0+fabs(g));//计算精度

　            p=g;

　            m=m+1;

　           h=(b-a)/m;//分割区间

　    }   

　    return g;

　}

　main()

　{

　    double a,b,s;

　    cout<<“请输入积分下限:“;

　    cin>>a;

　    cout<<“请输入积分上限:“;

　    cin>>b;

　    cout<<“㏑的真值为:“<<endl;

　    cout<<“1.098612289“<<endl;

　    /*龙贝格求积*/

　    s=Romberg( a, b, function);

　    cout<<“龙贝格求积公式:“<<endl;

　    cout<<setiosflags(ios::fixed)<<setprecision(14)<<s<<endl;

　     /*三点求积公式*/

　    s=ThreePointGaussLegendre( a, b, function);

　    cout<<“三点求积公式:“<<endl;

　    cout<<setiosflags(ios::fixed)<<setprecision(14)<<s<<endl;

　     /*五点求积公式*/

　    s=FivePointGaussLegendre( a, b, function);

　    cout<<“五点求积公式“<<endl;

　    cout<<setiosflags(ios::fixed)<<setprecision(14)<<s<<endl;

　     s=FivePointPrecisionGaussLegendre(a, b,function);

　    cout<<“控制精度五点求积公式“<<endl;

　    cout<<setiosflags(ios::fixed)<<setprecision(14)<<s<<endl;

　    return 0;

　} 2. 高斯－勒让德求积的程序

　2.1三点高斯勒让德公式的代码 function gl=f(str,a,b) x=zeros(3,1); y=zeros(3,1); x(1)=-sqrt(15)/5; x(2)=0; x(3)=sqrt(15)/5; for i=1:3

　 t=(b-a)/2*x(i)+(a+b)/2;

　 y(i)=eval(str);%exp(t)*sin(t);%此处为求积的函数,t为自变量 end gl=5/9*y(1)+8/9*y(2)+5/9*y(3);

　上面的代码保存为f.m文件，调用的时候如下

　f(＇t*2＇,-1,1)

　f(＇exp(t)*sin(t)＇,1,3)

　其中第一个参数为求积分的表达式，第二三个参数分别为

　积分的上下限。

　 2.2高斯-勒让德数值积分Matlab代码

　function [ql,Ak,xk]=guasslegendre(fun,a,b,n,tol) if nargin==1

　 a=-1;b=1;n=7;tol=1e-8; elseif nargin==3

　 n=7;tol=1e-8; elseif nargin==4

　 tol=1e-8; elseif nargin==2|nargin>5

　error(＇The Number of Input Arguments Is Wrong!＇); end syms x p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n)); tk=roots(p);

　Ak=zeros(n+1,1); for i=1:n+1

　 xkt=tk;

　 xkt(i)=[];

　 pn=poly(xkt);

　 fp=@(x)polyval(pn,x)/polyval(pn,tk(i));

　 Ak(i)=quadl(fp,-1,1,tol); % 求积系数 end xk=(b-a)/2*tk+(b+a)/2; fun=fcnchk(fun,＇vectorize＇); fx=fun(xk)*(b-a)/2; ql=sum(Ak.*fx);

　3.数值实验

　3.1

　用4点（n=3）的高斯——勒让德求积公式计算

　. 解：

　先将区间化为，由（1） .（1） 有 . 根据表4-7中n=3的节点及系数值可求得

　.

　（ 准确值 ）

　 3.2用的高斯-勒让德公式计算积分

　解：

　 令，则 用的高斯—勒让德公式计算积分

　用的高斯—勒让德公式计算积分

　3.2用四个节点的高斯―勒让德求积公式计算定积分，计算过程保留4位小数．

　 解 ：

　高斯－勒让德求积公式只求积分区间为[－1,1]上的积分问题．需作变换，令

　 ，当x=1时，u=1;当x=0时，u=－1．于是，

　＝

　 ＝

　3. 总结

　高斯―勒让德求积公式对定积分的计算拥有高精度的特点，但是这只存在于积分区间在[-1,1]上，区间的变大会导致精度的降低。因此，寻找精度更高，加速更快的算法是必要的。

　《参考文献》

　[1]《数值计算》

　张军、林瑛、钟竞辉 清华大学出版社

　2008 6 17

　[2]《数值分析》 陈晓江、黄樟灿· 科学出版社

　2010 7 10 [3]《数值分析原理》吴勃英

　科学出版社

　2009 7 23 [4] 复化两点Gauss-Legendre求积公式的外推算法 《桂林航天工业高等专科学校学报》2007年03期