【数值分析】高斯型求积公式，任意区间三点gauss求积公式，matlab实现

06-08 1564阅读

Gauss型求积公式

Gauss型求积公式定义

（图片来源网络，侵删）

∫ a b ρ ( x ) f ( x ) d x ≈ ∑ i = 1 n A i f ( x i ) \int_{ a }^{b} \rho(x)f(x) \mathrm dx \approx \sum_{i=1}^{ n}A_if(x_i) ∫abρ(x)f(x)dx≈i=1∑nAif(xi)

如果求积公式具有 2 n − 1 {2n-1} 2n−1 次代数精度，则称对应的节点 x 1 , x 2 , ⋯ , x n {x_1,x_2,\cdots ,x_n} x1,x2,⋯,xn 为Gauss点，此时求积公式称为Gauss型求积公式。

为了讨论方便，本节取 n {n} n 个节点，并记节点为 x 1 , x 2 , ⋯ x n {x_1,x_2,\cdots x_n} x1,x2,⋯xn ，从 1 {1} 1 开始取！同时，所讨论的积分均为带有权函数 ρ ( x ) {\rho(x)} ρ(x) 的积分。

插值型求积公式

∫ a b ρ ( x ) f ( x ) d x ≈ ∑ i = 1 n A i f ( x i ) \int_{ a }^{b} \rho(x)f(x) \mathrm dx \approx \sum_{i=1}^{ n}A_if(x_i) ∫abρ(x)f(x)dx≈i=1∑nAif(xi)

的代数精度最高不超过 2 n − 1 {2n-1} 2n−1 ，且达到 2 n − 1 {2n-1} 2n−1 时，所有求积系数为正。所以高斯公式是给定节点数下代数精度最高的求积公式。

构造Gauss型求积公式的步骤

对给定区间 [ a , b ] {[a,b]} [a,b] 及权函数 ρ ( x ) {\rho(x)} ρ(x) ，由Schmidt正交化过程构造正交多项式 P 0 ( x ) , P 1 ( x ) , ⋯ , P n ( x ) {P_0(x),P_1(x), \cdots ,P_n(x)} P0(x),P1(x),⋯,Pn(x)
求出 P n ( x ) {P_n(x)} Pn(x) 的 n {n} n 个零点 x 1 , x 2 , ⋯ , x n { x_1,x_2, \cdots, x_n} x1,x2,⋯,xn 即为Gauss点
计算求积系数 A i = ∫ a b ρ ( x ) l i ( x ) d x {A_i= \int_{ a }^{b} \rho(x)l_i(x) \mathrm dx } Ai=∫abρ(x)li(x)dx， i = 1 , 2 , ⋯ , n { i=1,2,\cdots,n } i=1,2,⋯,n， l i {l_i} li 为拉格朗日基函数

[!example]-

求计算积分 ∫ − 1 1 x 2 f ( x ) d x { \int_{ -1 }^{1} x^2f(x) \mathrm dx} ∫−11x2f(x)dx 的两点Gauss公式。

解：

ρ ( x ) = x 2 , n = 2 \rho(x)=x^2, n=2 ρ(x)=x2,n=2

首先按Schmidt正交化求出正交多项式：

P 0 ( x ) = 1 P 1 ( x ) = x − ( x , P 0 ( x ) ) ( P 0 ( x ) , P 0 ( x ) ) P 0 ( x ) = x − ∫ − 1 1 x 3 d x ∫ − 1 1 x 2 d x = x P 2 ( x ) = x 2 − ( x 2 , P 0 ( x ) ) ( P 0 ( x ) , P 0 ( x ) ) P 0 ( x ) − ( x 2 , P 1 ( x ) ) ( P 1 ( x ) , P 1 ( x ) ) P 1 ( x ) = x 2 − ∫ − 1 1 x 4 d x ∫ − 1 1 x 2 d x − ∫ − 1 1 x 5 d x ∫ − 1 1 x 4 d x x = x 2 − 3 5 \begin{align*} P_0(x)=&1 \\ \\ P_1(x)=&x- \frac{(x,P_0(x))}{(P_0(x),P_0(x))}P_0(x)=x- \frac{\int_{ -1 }^{1} x^3 \mathrm dx}{\int_{ -1 }^{1} x^2 \mathrm dx}=x \\ \\ P_2(x)=&x^2- \frac{(x^2,P_0(x))}{(P_0(x),P_0(x))}P_0(x)-\frac{(x^2,P_1(x))}{(P_1(x),P_1(x))}P_1(x) \\ \\ =&x^2- \frac{\int_{ -1 }^{1} x^4 \mathrm dx}{\int_{ -1 }^{1} x^2 \mathrm dx}- \frac{\int_{ -1 }^{1} x^5 \mathrm dx}{\int_{ -1 }^{1} x^4 \mathrm dx}x=x^2- \frac{3}{5} \end{align*} P0(x)=P1(x)=P2(x)==1x−(P0(x),P0(x))(x,P0(x))P0(x)=x−∫−11x2dx∫−11x3dx=xx2−(P0(x),P0(x))(x2,P0(x))P0(x)−(P1(x),P1(x))(x2,P1(x))P1(x)x2−∫−11x2dx∫−11x4dx−∫−11x4dx∫−11x5dxx=x2−53

再令 P 2 ( x ) = 0 {P_2(x)=0} P2(x)=0 求出Gauss点：

x 1 = − 3 5 , x 2 = 3 5 x_1=-\sqrt{\frac{3}{5}} \,\,,\,\, x_2=\sqrt{\frac{3}{5}} x1=−53 ,x2=53

最后计算求积系数：

A 1 = ∫ − 1 1 x 2 l 1 ( x ) d x = ∫ − 1 1 x 2 x − x 2 x 1 − x 2 d x = 1 3 A 2 = ∫ − 1 1 x 2 l 2 ( X ) d x = ∫ − 1 1 x 2 x − x 1 x 1 − x 2 d x = 1 3 \begin{align*} A_1=& \int_{ -1 }^{1} x^2l_1(x) \mathrm dx= \int_{ -1 }^{1} x^2 \frac{x-x_2}{x_1-x_2} \mathrm dx= \frac{1}{3}\\ \\ A_2=& \int_{ -1 }^{1} x^2l_2(X) \mathrm dx= \int_{ -1 }^{1} x^2 \frac{x-x_1}{x_1-x_2} \mathrm dx= \frac{1}{3} \end{align*} A1=A2=∫−11x2l1(x)dx=∫−11x2x1−x2x−x2dx=31∫−11x2l2(X)dx=∫−11x2x1−x2x−x1dx=31

所以两点Gauss公式为

∫ − 1 1 x 2 f ( x ) d x ≈ 1 3 [ f ( − 3 5 ) + f ( 3 5 ) ] \int_{ -1 }^{1} x^2f(x) \mathrm dx \approx \frac{1}{3}[f(- \sqrt{\frac{3}{5}})+f(\sqrt{\frac{3}{5}})] ∫−11x2f(x)dx≈31[f(−53 )+f(53 )]

4.1 Gauss-legendre勒让德求积公式

区间 [ − 1 , 1 ] {[-1,1]} [−1,1] ，权函数 ρ ( x ) = 1 {\rho(x)=1} ρ(x)=1 。

Gauss求积公式中在区间 [ − 1 , 1 ] {[-1,1]} [−1,1] 上权函数为 1 {1} 1 的Gauss点其实是确定的，所以可以通过把待求区间伸缩变换到 [ − 1 , 1 ] {[-1,1]} [−1,1] 间，再通过Gauss求积公式求解。

一点Gauss求积公式

∫ − 1 1 f ( x ) d x ≈ 2 f ( 0 ) \int_{ -1 }^{1} f(x) \mathrm dx \approx 2f(0) ∫−11f(x)dx≈2f(0)

两点Gauss求积公式

∫ − 1 1 f ( x ) d x ≈ f ( − 1 3 ) + f ( 1 3 ) \int_{ -1 }^{1} f(x) \mathrm dx \approx f(- \frac{1}{\sqrt{3}})+f(\frac{1}{\sqrt{3}}) ∫−11f(x)dx≈f(−3 1)+f(3 1)

三点Gauss求积公式

∫ − 1 1 f ( x ) d x ≈ 5 9 f ( − 3 5 ) + 8 9 f ( 0 ) + 5 9 f ( 3 5 ) \int_{ -1 }^{1} f(x) \mathrm dx \approx \frac{5}{9} f(- \sqrt{\frac{3}{5}})+ \frac{8}{9}f(0) +\frac{5}{9}f(\sqrt{\frac{3}{5}}) ∫−11f(x)dx≈95f(−53 )+98f(0)+95f(53 )

任意区间的三点高斯求积公式：

∫ a b f ( x ) d x ≈ ( b − a ) 2 ( 5 9 f ( ( a + b ) 2 − 3 5 ( b − a ) 2 ) + 8 9 f ( ( a + b ) 2 ) + 5 9 f ( ( a + b ) 2 + 3 5 ( b − a ) 2 ) ) \int_{ a }^{b} f(x) \mathrm dx \approx \frac{(b-a)}{2} \bigg( \frac{5}{9} f(\frac{(a+b)}{2} - \sqrt{\frac{3}{5}}\frac{(b-a)}{2})+ \frac{8}{9}f(\frac{(a+b)}{2}) +\frac{5}{9}f(\frac{(a+b)}{2} + \sqrt{\frac{3}{5}}\frac{(b-a)}{2}) \bigg) ∫abf(x)dx≈2(b−a)(95f(2(a+b)−53 2(b−a))+98f(2(a+b))+95f(2(a+b)+53 2(b−a)))

matlab实现

%% 三点高斯勒让德求积公式
% 输入函数，积分上界，积分下界
function I = gaussL3P(f,a,b)
    x = (a+b)/2+(b-a)/2*[-sqrt(0.6) 0 sqrt(0.6)];
    I = (b-a)/2*f(x)*[5 8 5]'/9;
end

四点Gauss求积公式

高斯点 : ± 0.8611363 , 高斯系数 : 0.3478548 高斯点:\pm0.8611363 \,\,,\,\, 高斯系数:0.3478548 高斯点:±0.8611363,高斯系数:0.3478548

高斯点 : ± 0.3399810 , 高斯系数 : 0.6521452 高斯点:\pm0.3399810 \,\,,\,\, 高斯系数:0.6521452 高斯点:±0.3399810,高斯系数:0.6521452

[!example]-

已知三点Gauss公式

∫ − 1 1 f ( x ) d x ≈ 5 9 f ( 0.6 ) + 8 9 f ( 0 ) + 5 9 f ( − 0.6 ) \int_{ -1 }^{1} f(x) \mathrm dx \approx \frac{5}{9}f(\sqrt{0.6})+ \frac{8}{9}f(0)+ \frac{5}{9}f(- \sqrt{0.6}) ∫−11f(x)dx≈95f(0.6 )+98f(0)+95f(−0.6 )

试用如上公式计算 ∫ 0.5 1 x d x { \int_{ 0.5 }^{1} \sqrt{x} \mathrm dx} ∫0.51x dx 的值。

解：换限，设

t = a x + b t=ax+b t=ax+b

{ 0.5 a + b = − 1 a + b = 1 ⇒ { a = 4 b = − 3 \begin{cases} 0.5a+b=-1 \\ \\ a+b=1 \end{cases} \Rightarrow \begin{cases} a=4 \\ \\ b=-3 \end{cases} ⎩ ⎨ ⎧0.5a+b=−1a+b=1⇒⎩ ⎨ ⎧a=4b=−3

∫ 0.5 1 f ( x ) d x = 1 4 ∫ − 1 1 f ( t + 3 4 ) d t ≈ 1 4 [ 5 9 f ( 0.6 + 3 4 ) + 8 9 f ( 0 + 3 4 ) + 5 9 f ( − 0.6 + 3 4 ) ] = 1 4 [ 5 9 0.6 + 3 4 + 8 9 0 + 3 4 + 5 9 − 0.6 + 3 4 ] = 0.4310 \begin{align*} \int_{ 0.5 }^{1} f(x) \mathrm dx=& \frac{1}{4} \int_{ -1 }^{1} f(\frac{t+3}{4}) \mathrm dt \\ \\ \approx& \frac{1}{4}[\frac{5}{9} f(\frac{\sqrt{0.6}+3}{4})+ \frac{8}{9}f(\frac{0+3}{4}) + \frac{5}{9}f(\frac{- \sqrt{0.6}+3}{4}) ] \\ \\ =&\frac{1}{4}[\frac{5}{9} \sqrt{\frac{\sqrt{0.6}+3}{4}}+ \frac{8}{9}\sqrt{\frac{0+3}{4}} + \frac{5}{9}\sqrt{\frac{- \sqrt{0.6}+3}{4}} ] \\ \\ =&0.4310 \end{align*} ∫0.51f(x)dx=≈==41∫−11f(4t+3)dt41[95f(40.6 +3)+98f(40+3)+95f(4−0.6 +3)]41[9540.6 +3 +9840+3 +954−0.6 +3 ]0.4310

4.2 变形的Gauss型求积公式

Gauss-Laguerre拉盖尔求积公式

区间 [ 0 , + ∞ ) {[0,+\infty)} [0,+∞) ，权函数 ρ ( x ) = e − x {\rho(x)=e^{-x}} ρ(x)=e−x

Gauss-Hermite求积公式

两点Gauss-Hermite求积公式

∫ − ∞ + ∞ e − x 2 f ( x ) d x ≈ π 2 f ( − 2 2 ) + π 2 f ( 2 2 ) \int_{ -\infty }^{+\infty} e^{-x^2}f(x) \mathrm dx \approx \frac{\sqrt{\pi}}{2} f(-\frac{\sqrt{2}}{2})+\frac{\sqrt{\pi}}{2} f(\frac{\sqrt{2}}{2}) ∫−∞+∞e−x2f(x)dx≈2π f(−22 )+2π f(22 )

三点Gauss-Hermite求积公式

∫ − ∞ + ∞ e − x 2 f ( x ) d x ≈ π 6 f ( − 6 2 ) + 2 π 3 f ( 0 ) + π 6 f ( 6 2 ) \int_{ -\infty }^{+\infty} e^{-x^2}f(x) \mathrm dx \approx \frac{\sqrt{\pi}}{6} f(-\frac{\sqrt{6}}{2})+ \frac{2 \sqrt{\pi}}{3} f(0)+\frac{\sqrt{\pi}}{6} f(\frac{\sqrt{6}}{2}) ∫−∞+∞e−x2f(x)dx≈6π f(−26 )+32π f(0)+6π f(26 )

区间 ( − ∞ , + ∞ ) {(-\infty,+\infty)} (−∞,+∞) ，权函数 ρ ( x ) = e − x 2 {\rho(x)=e^{-x^2}} ρ(x)=e−x2

VPS购买请点击我