Gauss型求积公式及其Matlab程序

为什么要引入数值积分

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

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

I ( f ) = ∫ a b f ( x ) d x I(f)=\int_{a}^{b} f(x) dx I(f)=abf(x)dx
的计算问题,从数学分析中知道,当已知 f ( x ) f(x) f(x)的原函数为 F ( x ) F(x) F(x)时,由牛顿-莱布尼兹公式,有
I ( f ) = ∫ a b f ( x ) d x = 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型求积公式即:

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

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

在这里插入图片描述
打代码太南了,我手写偷下懒。

Gauss公式的应用

Gauss公式主要应用为二点Gauss公式和三点Gauss公式
二点Gauss公式:当n=2时,有 x 1 = − 1 3 , x 2 = 1 3 , A 1 = A 2 = 1 x_1=-\sqrt\frac{1}{3},x_2=\sqrt\frac{1}{3},A_1=A_2=1 x1=31 ,x2=31 ,A1=A2=1,求积公式为 G 2 ( f ) = f ( − 1 3 ) + f ( 1 3 ) G_2(f)=f(-\sqrt\frac{1}{3})+f(\sqrt\frac{1}{3}) G2(f)=f(31 )+f(31 )
三点Gauss求积公式:当n=3时,有 x 1 = − 3 5 , x 2 = 0 , x 3 = 3 5 , A 1 = A 3 = 5 9 , A 2 = 8 9 x_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=0x3=53 ,A1=A3=95,A2=98,求积公式为 G 3 ( f ) = 5 9 f ( − 3 3 ) + 8 9 f ( 0 ) + 5 9 f ( 3 5 ) 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]上的定积分 ∫ a b f ( x ) d x \int_{a}^{b} f(x) dx abf(x)dx可以通过自变量的变换 x = 1 2 ( a + b + t ∗ ( b − a ) ) x=\frac{1}{2}(a+b+t*(b-a)) x=21(a+b+t(ba))将其转换为[-1,1]上的积分,然后使用Gauss积分公式计算其近似值。
例:
用两点Gauss公式计算 ∫ 0 1 s i n x x d x \int_{0}^{1}\frac{sinx}{x}dx 01xsinxdx
解:变换令x=0.5(t+1)
∫ 0 1 s i n x x d x = ∫ − 1 1 s i n 0.5 ( t + 1 ) t + 1 d x \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
t 0 = − 1 3 , t 1 = 1 3 t_0=\frac{-1}{\sqrt 3},t_1=\frac{1}{\sqrt 3} t0=3 1,t1=3 1
∫ 0 1 s i n x x d x = 0.5 [ s i n 0.5 ( t 0 + 1 ) 0.5 ( t 0 + 1 + s i n 0.5 ( t 1 + 1 ) 0.5 ( t 1 + 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 == 2
    w = [1,1];
    p = [-1/sqrt(3),1/sqrt(3)];
elseif n == 4
     w = [0.3478548451,0.3478548451,0.6521451549,0.6521451549];
     p = [0.8611363116,-0.8611363116,0.3399810436,-0.3399810436];
elseif n == 8
     w = [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));
  • 15
    点赞
  • 96
    收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:1024 设计师:我叫白小胖 返回首页
评论 2

打赏作者

凡尘红梦

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值