高斯—勒让德积分公式
摘要:
高斯—勒让德积分公式可以用较少节点数得到高精度的计算结果, 是现在现实生活中经常运用到的数值积分法。然而,当积分区间较大时,积分精度并不理想。
T he 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)
中,若取权函数
我们知道勒让
多项式,因此,勒让德多项式
德多项式是区间
上的正交
,区间为
,则得公式
的零点就是求积公式(上式) 的高斯点.形如
(上式) 的高斯公式特别地称为高斯-勒让德求积公式.
若取
令它对
准确成立,即可定出
.这样构造出的一点高斯-勒让
的零点
做节点构造求积公式
德求积公式是中矩形公式.再取式
令它对
都准确成立,有
的两个零点构造求积公
.
由此解出,从而得到两点高斯-勒让德求积公式
.
三点高斯-勒让德求积公式的形式是
.
如表列出高斯-勒让德求积公式的节点和系数.
公式(4.5.9)的余项由(4.5.8)得
,
这里是最高项系数为1的勒让德多项式,由(3.2.6)及(3.2.7)得
.
当时,有
.
它比辛普森公式余项还小,且比辛普森公式少算一个函数值.
当积分区间不是[-1,1]
,而是一般的区间时,只要做变换
可将
化为[-1,1],这时
.
对等式右端的积分即可使用高斯-勒让德求积公式.
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 #include #include
#define Precision1 0.[1**********]1 # 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
for(k=0;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
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
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
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
{
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>a;
cout>b;
cout
s=Romberg( a, b, function); cout
cout
/*三点求积公式*/
s=ThreePointGaussLegendre( a, b, function);
cout
cout
/*五点求积公式*/
s=FivePointGaussLegendre( a, b, function);
cout
cout
s=FivePointPrecisionGaussLegendre(a, b,function);
cout
cout
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)
其中第一个参数为求积分的表达式,第二三个参数分别为
11
积分的上下限。
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)的高斯——勒让德求积公式计算 π
I=⎰2
0x cos xdx . 2
解: 先将区间[0, π
2]化为[-1, 1],由(1)
12
⎰
I=b a f (x ) dx =b -a 2⎰1-1f (b -a 2+a +b 2) dt . (1) 有 ⎰1
-1(π4) (1+t ) cos 32π
4(1+t ) dt .
根据表4-7中n=3的节点及系数值可求得
3
I≈
∑A
k =0k f (x k ) ≈0. 467402.
( 准确值I=0. 467401 )
3.2用n =2, 3的高斯-勒让德公式计算积分
⎰
解:
I =31e sin xdx . x ⎰3
1e sin xdx . x
x ∈[1,3],令t =x -2,则t ∈[-1,1]
用n =2的高斯—勒让德公式计算积分
I ≈0.5555556⨯[f (-0.7745967) +f (0.7745967)]+0.8888889⨯f (0)
≈10.9484
用n =3的高斯—勒让德公式计算积分
I ≈0.3478548⨯[f (-0.8611363) +f (0.8611363)]
+0.6521452⨯[f (-0.3399810) +f (0.3399810)]
≈10.95014
3.2用四个节点的高斯―勒让德求积公式计算定积分⎰1
0+x d x ,计算过程保留4位小数.
解 :
高斯-勒让德求积公式只求积分区间为[-1,1]上的积分问题.需作变换,令
13
x =u
2+1
2,当x=1时,u=1;当x=0时,u=-1.于是,
1
⎰1
0+x d x =⎰2132-1+u 2d u
3
2-0. 8611
2+3
2+0. 8611
2) 1
=2[0. 3479⨯(
)]+0. 6521⨯(
1
232-0. 34002+32+0. 34002 =[0. 3479⨯2. 4236+0. 6521⨯2. 4455]=1. 2189
3. 总结
高斯―勒让德求积公式对定积分的计算拥有高精度的特点,但是这只存在于积分区间在[-1,1]上,区间的变大会导致精度的降低。因此,寻找精度更高,加速更快的算法是必要的。
《参考文献》
[1]《数值计算》 张军、林瑛、钟竞辉 清华大学出版社 2008 6 17
[2]《数值分析》 陈晓江、黄樟灿· 科学出版社 2010 7 10
[3]《数值分析原理》吴勃英 科学出版社 2009 7 23
[4] 复化两点Gauss-Legendre 求积公式的外推算法 《桂林航天工业高等专科学校学报》2007年03期
14
高斯—勒让德积分公式
摘要:
高斯—勒让德积分公式可以用较少节点数得到高精度的计算结果, 是现在现实生活中经常运用到的数值积分法。然而,当积分区间较大时,积分精度并不理想。
T he 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)
中,若取权函数
我们知道勒让
多项式,因此,勒让德多项式
德多项式是区间
上的正交
,区间为
,则得公式
的零点就是求积公式(上式) 的高斯点.形如
(上式) 的高斯公式特别地称为高斯-勒让德求积公式.
若取
令它对
准确成立,即可定出
.这样构造出的一点高斯-勒让
的零点
做节点构造求积公式
德求积公式是中矩形公式.再取式
令它对
都准确成立,有
的两个零点构造求积公
.
由此解出,从而得到两点高斯-勒让德求积公式
.
三点高斯-勒让德求积公式的形式是
.
如表列出高斯-勒让德求积公式的节点和系数.
公式(4.5.9)的余项由(4.5.8)得
,
这里是最高项系数为1的勒让德多项式,由(3.2.6)及(3.2.7)得
.
当时,有
.
它比辛普森公式余项还小,且比辛普森公式少算一个函数值.
当积分区间不是[-1,1]
,而是一般的区间时,只要做变换
可将
化为[-1,1],这时
.
对等式右端的积分即可使用高斯-勒让德求积公式.
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 #include #include
#define Precision1 0.[1**********]1 # 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
for(k=0;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
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
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
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
{
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>a;
cout>b;
cout
s=Romberg( a, b, function); cout
cout
/*三点求积公式*/
s=ThreePointGaussLegendre( a, b, function);
cout
cout
/*五点求积公式*/
s=FivePointGaussLegendre( a, b, function);
cout
cout
s=FivePointPrecisionGaussLegendre(a, b,function);
cout
cout
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)
其中第一个参数为求积分的表达式,第二三个参数分别为
11
积分的上下限。
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)的高斯——勒让德求积公式计算 π
I=⎰2
0x cos xdx . 2
解: 先将区间[0, π
2]化为[-1, 1],由(1)
12
⎰
I=b a f (x ) dx =b -a 2⎰1-1f (b -a 2+a +b 2) dt . (1) 有 ⎰1
-1(π4) (1+t ) cos 32π
4(1+t ) dt .
根据表4-7中n=3的节点及系数值可求得
3
I≈
∑A
k =0k f (x k ) ≈0. 467402.
( 准确值I=0. 467401 )
3.2用n =2, 3的高斯-勒让德公式计算积分
⎰
解:
I =31e sin xdx . x ⎰3
1e sin xdx . x
x ∈[1,3],令t =x -2,则t ∈[-1,1]
用n =2的高斯—勒让德公式计算积分
I ≈0.5555556⨯[f (-0.7745967) +f (0.7745967)]+0.8888889⨯f (0)
≈10.9484
用n =3的高斯—勒让德公式计算积分
I ≈0.3478548⨯[f (-0.8611363) +f (0.8611363)]
+0.6521452⨯[f (-0.3399810) +f (0.3399810)]
≈10.95014
3.2用四个节点的高斯―勒让德求积公式计算定积分⎰1
0+x d x ,计算过程保留4位小数.
解 :
高斯-勒让德求积公式只求积分区间为[-1,1]上的积分问题.需作变换,令
13
x =u
2+1
2,当x=1时,u=1;当x=0时,u=-1.于是,
1
⎰1
0+x d x =⎰2132-1+u 2d u
3
2-0. 8611
2+3
2+0. 8611
2) 1
=2[0. 3479⨯(
)]+0. 6521⨯(
1
232-0. 34002+32+0. 34002 =[0. 3479⨯2. 4236+0. 6521⨯2. 4455]=1. 2189
3. 总结
高斯―勒让德求积公式对定积分的计算拥有高精度的特点,但是这只存在于积分区间在[-1,1]上,区间的变大会导致精度的降低。因此,寻找精度更高,加速更快的算法是必要的。
《参考文献》
[1]《数值计算》 张军、林瑛、钟竞辉 清华大学出版社 2008 6 17
[2]《数值分析》 陈晓江、黄樟灿· 科学出版社 2010 7 10
[3]《数值分析原理》吴勃英 科学出版社 2009 7 23
[4] 复化两点Gauss-Legendre 求积公式的外推算法 《桂林航天工业高等专科学校学报》2007年03期
14