Gauss型求积公式及其Matlab程序

为什么要引入数值积分

由于Gauss型求积公式属于数值积分的内容，学东西总要知道它的来龙去脉，下面我简单介绍一下为什么要引入数值积分

给定函数f(x)∈C[a,b]f(x)\in C[a,b]f(x)∈C[a,b],考虑积分

I(f)=∫abf(x)dxI(f)=\int_{a}^{b} f(x) dxI(f)=∫abf(x)dx
的计算问题，从数学分析中知道，当已知f(x)f(x)f(x)的原函数为F(x)F(x)F(x)时，由牛顿-莱布尼兹公式，有
I(f)=∫abf(x)dx=F(b)−F(a)I(f)=\int_{a}^{b}f(x)dx=F(b)-F(a)I(f)=∫abf(x)dx=F(b)−F(a)
然而，在实际计算中，被积函数的f(x)f(x)f(x)的原函数经常无法用初等函数表示，但过于复杂。还有时，f(x)f(x)f(x)只在一些离散点上给出。在这样的情况下，就有必要借助数值方法来求I(f)I(f)I(f)的近似值。

Gauss型求积公式的定义

从正交多项式理论可知，在区间[a,b]上，对给定的权函数ρ(x)\rho(x)ρ(x),存在正交多项式系ωk(x)k=0∞{\omega_{k}(x)}_{k=0}^{\infty}ωk(x)k=0∞并且可以将其构造出来。进一步，已经证明了ωk(x)\omega_{k}(x)ωk(x)在[a,b]上恰有k个相异的根，取ωn+1(x)\omega_{n+1}(x)ωn+1(x)的n+1个根为求积结点，构造插值型求积公式，将得到具有2n+1次代数精度的求积公式，称如此构造出来的求积公式为Gauss型求积公式即：

如果求积结点x0,x1,....,xnx_{0},x_{1},....,x_{n}x0,x1,....,xn,使插值型求积公式∫−11f(x)dx≈∑k=0nAkf(xk)\int_{-1}^{1}f(x)dx\approx\sum_{k=0}^{n}A_kf(x_k)∫−11f(x)dx≈∑k=0nAkf(xk)的代数精度为2n+1，则称该求积公式为Gauss型求积公式，称这些求积结点为Gauss点。

特别地，Gauss求积公式当[a,b]=[-1,1]时，取ρ(x)\rho(x)ρ(x)=1,Gauss型求积公式称为Gauss求积公式，此时的求积结点
xk(k=1,2,3....)x_k(k=1,2,3....)xk(k=1,2,3....)为n次Legendre多项式的根。
且Ak=∫−11ω(x)(x−xk)ω′(xk)dxA_k=\int_{-1}^{1} \frac{\omega(x)}{(x-x_k)\omega'(x_k)}dxAk=∫−11(x−xk)ω′(xk)ω(x)dx
而Legendre多项式为：

打代码太南了，我手写偷下懒。

Gauss公式的应用

Gauss公式主要应用为二点Gauss公式和三点Gauss公式
二点Gauss公式：当n=2时，有 x1=−13,x2=13,A1=A2=1x_1=-\sqrt\frac{1}{3},x_2=\sqrt\frac{1}{3},A_1=A_2=1x1=−31,x2=31,A1=A2=1,求积公式为G2(f)=f(−13)+f(13)G_2(f)=f(-\sqrt\frac{1}{3})+f(\sqrt\frac{1}{3})G2(f)=f(−31)+f(31)
三点Gauss求积公式：当n=3时，有x1=−35,x2=0，x3=35,A1=A3=59,A2=89x_1=-\sqrt\frac{3}{5},x_2=0，x_3=\sqrt\frac{3}{5},A_1=A_3=\frac{5}{9},A_2=\frac{8}{9}x1=−53,x2=0，x3=53,A1=A3=95,A2=98,求积公式为G3(f)=59f(−33)+89f(0)+59f(35)G_3(f)=\frac{5}{9}f(-\sqrt\frac{3}{3})+\frac{8}{9}f(0)+\frac{5}{9}f(\sqrt\frac{3}{5})G3(f)=95f(−33)+98f(0)+95f(53)
为了计算[a,b]上的定积分∫abf(x)dx\int_{a}^{b} f(x) dx∫abf(x)dx可以通过自变量的变换x=12(a+b+t∗(b−a))x=\frac{1}{2}(a+b+t*(b-a))x=21(a+b+t∗(b−a))将其转换为[-1,1]上的积分，然后使用Gauss积分公式计算其近似值。
例：
用两点Gauss公式计算∫01sinxxdx\int_{0}^{1}\frac{sinx}{x}dx∫01xsinxdx
解：变换令x=0.5(t+1)
∫01sinxxdx=∫−11sin0.5(t+1)t+1dx\int_{0}^{1}\frac{sinx}{x}dx=\int_{-1}^{1}\frac{sin0.5(t+1)}{t+1}dx∫01xsinxdx=∫−11t+1sin0.5(t+1)dx
取t0=−13,t1=13t_0=\frac{-1}{\sqrt 3},t_1=\frac{1}{\sqrt 3}t0=3−1,t1=31
∫01sinxxdx=0.5[sin0.5(t0+1)0.5(t0+1+sin0.5(t1+1)0.5(t1+1]\int_{0}^{1}\frac{sinx}{x}dx=0.5[\frac{sin0.5(t_0+1)}{0.5(t_0+1}+\frac{sin0.5(t_1+1)}{0.5(t_1+1}]∫01xsinxdx=0.5[0.5(t0+1sin0.5(t0+1)+0.5(t1+1sin0.5(t1+1)]

Matlab程序

function [w,p] = Gauss_point_1D(n,a,b)
% Gauss quarature point on [-1,1]
if n == 2w = [1,1];p = [-1/sqrt(3),1/sqrt(3)];
elseif n == 4w = [0.3478548451,0.3478548451,0.6521451549,0.6521451549];p = [0.8611363116,-0.8611363116,0.3399810436,-0.3399810436];
elseif n == 8w = [0.1012285363,0.1012285363,0.2223810345,0.2223810345,0.3137066459,0.3137066459,0.3626837834,0.3626837834];p = [0.9602898565,-0.9602898565,0.7966664774,-0.7966664774,0.5255324099,-0.5255324099,0.1834346425,-0.1834346425];
end% Gauss quarature point on [a,b]
w = 0.5*(b-a)*w;
p = 0.5*(b-a)*p+0.5*(b+a);function q = Gauss_quadrature(fun,a,b,n)
% Gauss quadrature
[w,p] = Gauss_point_1D(n,a,b);
q = sum(w.*fun(p));