【数值计算之二】数值积分之牛顿——科斯特公式：梯形、辛普森、辛普森3/8和布尔 高斯积分公式：勒让德、切比雪夫、拉盖尔和埃尔米特
import numpy as np
from scipy.integrate import quad
from sympy import *
init_printing()
import matplotlib.pyplot as plt
数值积分
数值积分的本质都是采用了微积分的思想：化整为零（将积分区间离散化），以直代曲（基于积分中值定理，用函数值的加权平均数近似表示区间的平均高度）
- 闭型牛顿——科斯特公式：梯形公式、辛普森(simpson)公式和布尔(Boole)公式，都是选择等距节点；
- 高斯积分公式：高斯——勒让德、高斯——切比雪夫、高斯——拉盖尔和高斯——埃尔米特。高斯——勒让德(Gauss-Legendre)公式选择某些勒让德多项式的根作为不等距节点。类似的还有高斯——切比雪夫积分、高斯——拉盖尔积分和高斯——埃尔米特积分等等。
- 积分中值定理——平均高度的各种近似法
共同点：
(a) 都依赖于步长（梯形宽度），权重一般呈中间高，两边低的对称分布。
(b) 用基于各个节点x0,x1,...,xMx_0, x_1,...,x_Mx0,x1,...,xM的f(x)f(x)f(x)的拉格朗日逼近多项式PM(x)P_M(x)PM(x)来代替被积函数f(x)f(x)f(x)，本质上就是取各个节点的"加权平均数"，这里的权重之和不一定为1。不同点
(a) 闭型牛顿——科斯特公式选择等距节点，这限制了求积公式的代数精度；而高斯——勒让德(Gauss-Legendre)公式则取消了这个限制条件，因此能够极大地提高了代数精度。从数值分析可以知道, 高斯方法的确比梯形法, 辛普森, 龙格库塔法有优越性: 对被积函数的拟合精度高, 收敛快。
1.闭型牛顿——科斯特公式
假设积分区间固定为[a,b][a,b][a,b]，那么不同公式要使用不一样的步长hhh。
1.梯形公式：2个点，将积分区间分为111等分,精度为 n=1n=1n=1，步长 h=b−ah=b-ah=b−a
2.辛普森公式：3个点,将积分区间分为222等分,精度为 n=3n=3n=3,步长 h=b−a2h=\frac{b-a}{2}h=2b−a
3.辛普森 38\frac{3}{8}83 公式：4个点,将积分区间分为333等分,精度为 n=3n=3n=3,步长 h=b−a3h=\frac{b-a}{3}h=3b−a
4.布尔公式：5个点,将积分区间分为444等分,精度为 n=5n=5n=5, 步长 h=b−a4h=\frac{b-a}{4}h=4b−a
5.组合梯形：将积分区间分为 MMM 等分,MMM 为正整数，共(M+1)(M+1)(M+1)个点，步长为 h=b−aMh=\frac{b-a}{M}h=Mb−a，误差的阶为 O(h2)O(h^{2})O(h2)
6.组合辛普森公式：将积分区间分为2M2M2M等分，共(2M+1)(2M+1)(2M+1)个点，步长为 h=b−a2Mh=\frac{b-a}{2M}h=2Mb−a,误差的阶为 O(h4)O(h^{4})O(h4)
各公式的具体表达式如下：
Python 3.6 代码
def trapezoid(h ,f):Integral_appro = (h / 2) * (f[0] + f[1]) return Integral_approdef simpson(h, f):Integral_appro = (h / 3) * (f[0] + 4 * f[1] + f[2])return Integral_approdef simpson_3_8(h,f):Integral_appro = (3/8 * h) * (f[0] + 3 * f[1] + 3 * f[2] + f[3])return Integral_approdef boole(h, f):Integral_appro = (2/45 * h) * (7 * f[0] + 32 * f[1] + 12 * f[2] + 32 * f[3] + 7 * f[4])return Integral_appro闭型牛顿——科斯特公式假设积分区间为[0,1]，那么各公式的步长 h 依次为1，1/2， 1/3， 1/4, 也就是分别将区间等分为1、2、3、4份，取各个节点的函数值 f(xi)f(x_i)f(xi) 的"加权平均数" w(xi)f(xi)w(x_i)f(x_i)w(xi)f(xi) 作为整个区间平均高度f(e)f(e)f(e)的近似
def func(x):f = 1 + np.exp(-x) * np.sin(4 * x)return f
I_precise = quad(func,0,1)
I_precise
(1.3082506046426685,1.4524499432540732e−14)\left ( 1.3082506046426685, \quad 1.4524499432540732e-14\right )(1.3082506046426685,1.4524499432540732e−14)
x_s = Symbol('x_s', real=True)
f_s = 1 + exp(-x_s) * sin(4 * x_s)
f_s
1+e−xssin(4xs)1 + e^{- x_{s}} \sin{\left (4 x_{s} \right )}1+e−xssin(4xs)
Inte_f_symbol = Integral(f_s,(x_s,0, 1)) # 积分表达式（只是符号，没有求出原函数）
Inte_f_s = integrate(f_s, x_s) # 积分表达式（求出原函数）Inte_f_s_n1 = integrate(f_s, (x_s,0,1)) # 定积分
Inte_f_s_n2 = integrate(f_s, (x_s,0,1)).evalf() # 定积分，并将结果化成小数Inte_f_symbol,Inte_f_s,Inte_f_s_n1,Inte_f_s_n2
(∫01(1+e−xssin(4xs)) dxs,xs−e−xs17sin(4xs)−417e−xscos(4xs),−sin(4)17e−417ecos(4)+2117,1.30825060464267)\left ( \int_{0}^{1} \left(1 + e^{- x_{s}} \sin{\left (4 x_{s} \right )}\right)\, dx_{s}, \quad x_{s} - \frac{e^{- x_{s}}}{17} \sin{\left (4 x_{s} \right )} - \frac{4}{17} e^{- x_{s}} \cos{\left (4 x_{s} \right )}, \quad - \frac{\sin{\left (4 \right )}}{17 e} - \frac{4}{17 e} \cos{\left (4 \right )} + \frac{21}{17}, \quad 1.30825060464267\right )(∫01(1+e−xssin(4xs))dxs,xs−17e−xssin(4xs)−174e−xscos(4xs),−17esin(4)−17e4cos(4)+1721,1.30825060464267)
# 假设积分区间为[0,1]，那么各公式的步长 h 依次为1，1/2， 1/3， 1/4,分别将区间分为1、2、3、4份，
# 取各个节点的函数值的加权平均数作为整个区间平均高度的近似
I_trapezoid = trapezoid(1,[func(0), func(1)])
I_simpson = simpson(1/2, [func(0), func(1/2), func(1)])
I_simpson_3_8 = simpson_3_8(1/3, [func(0), func(1/3), func(2/3), func(1)])
I_boole = boole(1/4, [func(0), func(1/4), func(2/4), func(3/4),func(1)])
I_trapezoid, I_simpson, I_simpson_3_8, I_boole
(0.8607939604744832,1.3212758322698814,1.3143968149336276,1.3085919215646966)\left ( 0.8607939604744832, \quad 1.3212758322698814, \quad 1.3143968149336276, \quad 1.3085919215646966\right )(0.8607939604744832,1.3212758322698814,1.3143968149336276,1.3085919215646966)
def func(x):f = 2 + np.sin(2 * np.sqrt(x))return f# M 等分
def composite_trapezoid(f, a, b, M):if M < 1:print('M must be larger than or equal to 1')returnh = (b - a) / Ms1 = h / 2 * (f(a) + f(b))s2 = 0for k in range(1,M):x = a + h * ks2 += f(x)s = s1 + s2 * hreturn s# 2M 等分
def composite_simpson(f, a, b, M):if M < 1:print('M must be larger than or equal to 1')returnh = (b - a) / (2 * M)s1 = h / 3 * (f(a) + f(b))s2 = 0# k = 1, 2,...,Mfor k in range(1,M+1):x = a + h * (2 * k - 1)s2 += f(x) s2 = 4 * h / 3 * s2s3 = 0# k = 1,2,...,M-1for k in range(1, M):x = a + h * k * 2s3 += f(x)s3 = 2 * h / 3 * s3s = s1 + s2 + s3return sI_true = quad(func,1,6)
I_tra_M_10 = composite_trapezoid(func, 1, 6, 10)
I_tra_M_20 = composite_trapezoid(func, 1, 6, 20)
I_tra_M_40 = composite_trapezoid(func, 1, 6, 40)
I_tra_M_60 = composite_trapezoid(func, 1, 6, 60)
I_tra_M_80 = composite_trapezoid(func, 1, 6, 80)
I_tra_M_100 = composite_trapezoid(func, 1, 6, 100)
I_true,I_tra_M_10,I_tra_M_20,I_tra_M_40,I_tra_M_60,I_tra_M_80,I_tra_M_100
((8.183479207662728,1.79863015057564e−10),8.193854565172531,8.186049263770313,8.184120191790313,8.183763962635512,8.183639357318619,8.183581696027387)\left ( \left ( 8.183479207662728, \quad 1.79863015057564e-10\right ), \quad 8.193854565172531, \quad 8.186049263770313, \quad 8.184120191790313, \quad 8.183763962635512, \quad 8.183639357318619, \quad 8.183581696027387\right )((8.183479207662728,1.79863015057564e−10),8.193854565172531,8.186049263770313,8.184120191790313,8.183763962635512,8.183639357318619,8.183581696027387)
I_sim_M_10 = composite_simpson(func, 1, 6, 10)
I_sim_M_20 = composite_simpson(func, 1, 6, 20)
I_sim_M_10,I_sim_M_20
(8.183447496636239,8.183477167796982)\left ( 8.183447496636239, \quad 8.183477167796982\right )(8.183447496636239,8.183477167796982)
2. 高斯——勒让德公式
2.1 高斯积分
高斯积分法是精度最高的插值型数值积分，具有2n−12n-12n−1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到。
高斯积分有很多种，高斯——勒让德只是其中常用的一种。如果对具体的推导过程及其它类型的Gauss积分有兴趣，请查阅参考资料 [4] 高斯——勒让德公式–中南大学 和
[5] 各种类型的高斯积分 ，讲得很详细！
如果你想更进一步，请看这里！
目前最高效的数值积分方法：Gauss–Kronrod quadrature formula及其变种，这是一种自适应的数值积分方法，在高斯——勒让德的nnn个节点中，根据Stieltjes多项式增加了n+1n+1n+1个节点，因此一共有2n+12n+12n+1个节点，特别适合于高度振荡的函数。
详情请看[3] wiki:Gauss–Kronrod quadrature formula，目前最高效的实现为Laurie(1997)的论文[8] 及后来Calvetti 等人(2000)提出的改进版本[9]（在百度学术上都能找到可免费下载的资源）
Virginia Tech 的 John Burkardt 老师在2016年给出了 Gauss–Kronrod 积分的 Python 3实现：http://people.sc.fsu.edu/~jburkardt/py_src/kronrod/kronrod.html
下面摘录Calvetti 等人(2000)的文章中一些关键的定义：
好啦，来聊些简单的东西。
首先看看高斯积分的定义：
注意：这张PPT是我从参考资料中随手截取的，其实它取了n+1n+1n+1个点，因此精度为2n+12n+12n+1，而下面的讨论都是取nnn个点，别蒙圈了兄弟
注：关于正交多项式的定义，在一个国外的课件找到另外一种定义，暂时还没搞清楚哪一种更准确，先记录一下：
感叹一下：百度文库进步了，连文献和国外的PPT都有！
https://wenku.baidu.com/view/1d88b420bcd126fff7050b0c.html?rec_flag=default
- 不同的权函数 ρ(x)\rho(x)ρ(x)对应不同的正交多项式 P(x)P(x)P(x)。
- 2.1.1 高斯——勒让德积分
高斯——勒让德积分的权函数 ρ(x)=1\rho(x)=1ρ(x)=1，对应的正交多项式P(x)P(x)P(x)为勒让德多项式 Pn(x)P_n(x)Pn(x)，积分区间为[−1,1][-1,1][−1,1]。
高斯——勒让德积分的本质就是被积函数 f(x)f(x)f(x) 在区间 (−1,1)(-1,1)(−1,1) 内的积分近似于区间 (−1,1)(-1,1)(−1,1) 内的 nnn个勒让德——多项式零点的函数值 f(xk)f(x_k)f(xk) 的加权平均数 sum(wkf(xk))sum(w_k f(x_k))sum(wkf(xk))，权重之和不一定为1，权重依然呈现中间高，两边低的趋势。
- 2.1.2 高斯——切比雪夫积分
高斯——切比雪夫积分的权函数 ρ(x)=11−x2\rho(x)=\frac{1}{1-x^2}ρ(x)=1−x21，对应的正交多项式P(x)P(x)P(x)为切比雪夫多项式 Tn(x)T_n(x)Tn(x)，积分区间为[−1,1][-1,1][−1,1]。**
- 2.1.3 高斯——拉盖尔积分
高斯——拉盖尔积分的权函数 ρ(x)=e−x\rho(x)=e^{-x}ρ(x)=e−x，对应的正交多项式P(x)P(x)P(x)为拉盖尔多项式 Ln(x)L_n(x)Ln(x)，积分区间为[0,+∞][0,+\infty][0,+∞]。
- 2.1.4 高斯——埃尔米特积分
高斯——埃尔米特积分的权函数 ρ(x)=e−x2\rho(x)=e^{-x^{2}}ρ(x)=e−x2，对应的正交多项式P(x)P(x)P(x)为埃尔米特多项式 Hn(x)H_n(x)Hn(x),积分区间为 [−∞,+∞][-\infty,+\infty][−∞,+∞]。
numpy有对应的模块：
Gauss-Legendre:
节点和权值：numpy.polynomial.legendre.leggauss
权函数ρ(x)=1\rho(x)=1ρ(x)=1：numpy.polynomial.legendre.legweightGauss-Chebyshev:
节点和权值：numpy.polynomial.chebyshev.chebgauss
权函数ρ(x)=11−x2\rho(x)=\frac{1}{1-x^2}ρ(x)=1−x21：numpy.polynomial.chebyshev.chebweightGauss-Laguerre:
节点和权值：numpy.polynomial.laguerre.laggauss
权函数ρ(x)=e−x\rho(x)=e^{-x}ρ(x)=e−x：numpy.polynomial.laguerre.lagweightGauss-Hermite:
节点和权值：numpy.polynomial.hermite.hermgauss
权函数ρ(x)=e−x2\rho(x)=e^{-x^{2}}ρ(x)=e−x2：numpy.polynomial.hermite.hermweight
- 对于任意有限区间 [a,b],a,b≠±∞[a,b], a,b≠±\infty[a,b],a,b̸=±∞ 内的积分，可以借助线性变换将积分区间变换到[−1,1][-1,1][−1,1]，从而利用高斯——勒让德积分和高斯——切比雪夫积分进行计算。
2.2 高斯-勒让德积分
高斯-勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定，让函数f(x)f(x)f(x) 依次取i=0,1,2....ni=0,1,2....ni=0,1,2....n次多项式，使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n−12n−12n−1，而且是最高的。通常运用的是(−1,1)(−1,1)(−1,1)的积分节点和积分系数，其他积分域是通过一定的变换变换到-1到1之间积分。
勒让德多项式
首先复习一下勒让德多项式：
一般NNN点的高斯——勒让德公式对于次数小于等于2N−12N-12N−1次的多项式是精确的，即GN(f)G_N(f)GN(f)精度为n=2N−1n=2N-1n=2N−1,误差项E[f]E[f]E[f]包含2N2N2N阶导数。
Q：为什么NNN点的高斯——勒让德公式对于次数小于等于2N−12N-12N−1次的多项式是精确的呢？
A：因为NNN点的高斯——勒让德公式包含NNN个未知节点和相应的NNN个未知权值，一个有2N2N2N个未知数。我们可以从0阶多项式开始构造方程，则2N2N2N个未知数需要2N2N2N个等式方程，需要多项式 xc,c=0,1,2,...2N−1x^c, c=0,1,2,...2N-1xc,c=0,1,2,...2N−1，所以对次数小于等于2N−12N-12N−1次的多项式是可以取等号的，因此是精确的，而对2N−12N-12N−1次及以上的多项式则只能取约等于号，存在误差。
因此222点的高斯——勒让德多项式对于333次的多项式是精确的，333点的高斯——勒让德多项式对于555次的多项式是精确的。
因此，222点的高斯——勒让德多项式的误差项E[f]E[f]E[f]包含444阶导数，333点的高斯——勒让德多项式的误差项E[f]E[f]E[f]包含666阶导数。
NNN点的高斯——勒让德公式的节点就是nnn阶勒让德多项式Pn(x)P_n(x)Pn(x)的根。
想要在区间t∈[a,b]t∈[a,b]t∈[a,b]上使用高斯——勒让德公式，只需作变量替换t=(a+b)/2+(b−a)/2∗xt=(a+b)/2 + (b-a)/2 * xt=(a+b)/2+(b−a)/2∗x，节点和权值必须从高斯——勒让德节点与权值表中获得。
以下是Python3.6的代码实现。
2.1 用Sympy和numpy求勒让德多项式的节点
n阶勒让德多项式
x = Symbol('x', real=True)
n = Symbol('n', real=True)
legendre(n,x)
Pn(x)P_{n}\left(x\right)Pn(x)
0——4阶勒让德多项式及多项式的根
legendre(0,x), solve(legendre(0,x),x) # 0阶
(1,[])\left ( 1, \quad \left [ \right ]\right )(1,[])
legendre(1,x), solve(legendre(1,x),x) # 1阶
(x,[0])\left ( x, \quad \left [ 0\right ]\right )(x,[0])
legendre(2,x), solve(legendre(2,x),x) # 2阶
(3x22−12,[−33,33])\left ( \frac{3 x^{2}}{2} - \frac{1}{2}, \quad \left [ - \frac{\sqrt{3}}{3}, \quad \frac{\sqrt{3}}{3}\right ]\right )(23x2−21,[−33,33])
legendre(3,x), solve(legendre(3,x),x) # 3阶
(5x32−3x2,[0,−155,155])\left ( \frac{5 x^{3}}{2} - \frac{3 x}{2}, \quad \left [ 0, \quad - \frac{\sqrt{15}}{5}, \quad \frac{\sqrt{15}}{5}\right ]\right )(25x3−23x,[0,−515,515])
legendre(4,x),solve(legendre(4,x),x) # 4阶
(35x48−15x24+38,[−−23035+37,−23035+37,−23035+37,23035+37])\left ( \frac{35 x^{4}}{8} - \frac{15 x^{2}}{4} + \frac{3}{8}, \quad \left [ - \sqrt{- \frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad \sqrt{- \frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad - \sqrt{\frac{2 \sqrt{30}}{35} + \frac{3}{7}}, \quad \sqrt{\frac{2 \sqrt{30}}{35} + \frac{3}{7}}\right ]\right )⎝⎛835x4−415x2+83,⎣⎡−−35230+73,−35230+73,−35230+73,35230+73⎦⎤⎠⎞
求 f=1/xf = 1 / xf=1/x在[1,5]上的积分
def func(x):f = 1 / xreturn f# nodes是标准的勒让德多项式的根
# T是替换变量
# 对于不同阶数的勒让德多项式，nodes和weights是不相同的
def gauss_lengendre(fun, a,b, nodes,weights):nodes = np.array(nodes)weights = np.array(weights)T = (a+b)/2 + (b-a)/2 * nodes # 变量替换quad = (b-a)/2 * np.sum(weights * fun(T)) # element-wise multiplicationreturn quad
精确值
I_true = quad(func,1,5)
I_true
(1.6094379124341014,3.6599536780638603e−09)\left ( 1.6094379124341014, \quad 3.6599536780638603e-09\right )(1.6094379124341014,3.6599536780638603e−09)
用3点勒让德多项式近似
le_3 = legendre(3,x)
nodes_3 = solve(le_3,x)
nodes_3
[0,−155,155]\left [ 0, \quad - \frac{\sqrt{15}}{5}, \quad \frac{\sqrt{15}}{5}\right ][0,−515,515]
weights = np.array([5/9, 8/9, 5/9])
nodes = np.array([-np.sqrt(15)/5, 0, np.sqrt(15)/5])
a = 1
b = 5
I_n_3 = gauss_lengendre(func,a,b, nodes, weights)
error_n_3 = abs(I_true[0] - I_n_3)
I_n_3,error_n_3
(1.6026936026936027,0.006744309740498666)\left ( 1.6026936026936027, \quad 0.006744309740498666\right )(1.6026936026936027,0.006744309740498666)
用8点勒让德多项式近似
le_8 = legendre(8,x)
nodes_8 = solve(le_8,x) # 求出根节点
nodes_8_f = [i.evalf() for i in nodes_8] # list nodes_8_f 的每个元素都是 sympy.core.numbers.Float
nodes_8_arr = np.array(nodes_8_f, dtype=float) # 转换成ndarray，将数据类型改为float
nodes_8_arr
array([-0.18343464, 0.18343464, -0.52553241, 0.52553241, -0.79666648,0.79666648, -0.96028986, 0.96028986])
weights_8_arr = np.array([0.3626837834, 0.3626837834, 0.3137066459, 0.3137066459, 0.2223810345, 0.2223810345, 0.1012285363,0.1012285363], dtype=np.float64)
weights_8_arr
array([0.36268378, 0.36268378, 0.31370665, 0.31370665, 0.22238103,0.22238103, 0.10122854, 0.10122854])
a = 1
b = 5
I_n_8 = gauss_lengendre(func,a,b, nodes_8_arr, weights_8_arr)
error_n_8 = abs(I_true[0] - I_n_8)
I_n_8, error_n_8
(1.609437439052705,4.733813963042621e−07)\left ( 1.609437439052705, \quad 4.733813963042621e-07\right )(1.609437439052705,4.733813963042621e−07)
2.2 用numpy.polynomial.legendre.leggauss(n)一键求出n阶高斯——勒让德节点和权值
numpy.polynomial.legendre.leggauss(n) 函数可以一键求出nnn阶高斯——勒让德节点和权值，在n<=100n<=100n<=100的范围内，这个函数的结果是准确的。
n_100, w_100 = np.polynomial.legendre.leggauss(100)
n_100, w_100
(array([-0.99971373, -0.99849195, -0.99629513, -0.99312494, -0.9889844 ,-0.98387754, -0.97780936, -0.97078578, -0.96281365, -0.95390078,-0.94405587, -0.93328854, -0.9216093 , -0.90902957, -0.89556164,-0.88121868, -0.86601469, -0.84996453, -0.83308388, -0.81538924,-0.79689789, -0.77762791, -0.75759812, -0.73682809, -0.71533812,-0.6931492 , -0.67028302, -0.64676191, -0.62260886, -0.59784747,-0.57250193, -0.54659701, -0.52015802, -0.49321079, -0.46578165,-0.4378974 , -0.40958529, -0.38087298, -0.35178853, -0.32236034,-0.29261719, -0.26258812, -0.23230248, -0.20178986, -0.17108008,-0.14020314, -0.1091892 , -0.07806858, -0.04687168, -0.01562898,0.01562898, 0.04687168, 0.07806858, 0.1091892 , 0.14020314,0.17108008, 0.20178986, 0.23230248, 0.26258812, 0.29261719,0.32236034, 0.35178853, 0.38087298, 0.40958529, 0.4378974 ,0.46578165, 0.49321079, 0.52015802, 0.54659701, 0.57250193,0.59784747, 0.62260886, 0.64676191, 0.67028302, 0.6931492 ,0.71533812, 0.73682809, 0.75759812, 0.77762791, 0.79689789,0.81538924, 0.83308388, 0.84996453, 0.86601469, 0.88121868,0.89556164, 0.90902957, 0.9216093 , 0.93328854, 0.94405587,0.95390078, 0.96281365, 0.97078578, 0.97780936, 0.98387754,0.9889844 , 0.99312494, 0.99629513, 0.99849195, 0.99971373]),array([0.00073463, 0.00170939, 0.00268393, 0.00365596, 0.00462445,0.00558843, 0.00654695, 0.00749907, 0.00844387, 0.00938042,0.0103078 , 0.01122511, 0.01213146, 0.01302595, 0.01390771,0.01477588, 0.01562962, 0.01646809, 0.01729046, 0.01809594,0.01888374, 0.01965309, 0.02040323, 0.02113344, 0.021843 ,0.02253122, 0.02319742, 0.02384096, 0.0244612 , 0.02505754,0.0256294 , 0.02617622, 0.02669746, 0.02719261, 0.0276612 ,0.02810276, 0.02851685, 0.02890309, 0.02926108, 0.02959049,0.02989098, 0.03016227, 0.03040408, 0.03061619, 0.03079838,0.03095048, 0.03107234, 0.03116384, 0.03122488, 0.03125542,0.03125542, 0.03122488, 0.03116384, 0.03107234, 0.03095048,0.03079838, 0.03061619, 0.03040408, 0.03016227, 0.02989098,0.02959049, 0.02926108, 0.02890309, 0.02851685, 0.02810276,0.0276612 , 0.02719261, 0.02669746, 0.02617622, 0.0256294 ,0.02505754, 0.0244612 , 0.02384096, 0.02319742, 0.02253122,0.021843 , 0.02113344, 0.02040323, 0.01965309, 0.01888374,0.01809594, 0.01729046, 0.01646809, 0.01562962, 0.01477588,0.01390771, 0.01302595, 0.01213146, 0.01122511, 0.0103078 ,0.00938042, 0.00844387, 0.00749907, 0.00654695, 0.00558843,0.00462445, 0.00365596, 0.00268393, 0.00170939, 0.00073463]))
I_n_100 = gauss_lengendre(func,a,b, n_100, w_100)
error_n_100 = abs(I_true[0] - I_n_100)
I_n_100, error_n_100
(1.6094379124340965,4.884981308350689e−15)\left ( 1.6094379124340965, \quad 4.884981308350689e-15\right )(1.6094379124340965,4.884981308350689e−15)
error_n_3,error_n_8, error_n_100
(0.006744309740498666,4.733813963042621e−07,4.884981308350689e−15)\left ( 0.006744309740498666, \quad 4.733813963042621e-07, \quad 4.884981308350689e-15\right )(0.006744309740498666,4.733813963042621e−07,4.884981308350689e−15)
比较不同阶数的计算误差，可见高阶的一般能取得较高精度。
[1]高斯——勒让德积分中不同阶数下最大高斯节点间距的关系 节点数N与最大间距呈反相关;帖子给出了2-8阶勒让德多项式的权值
[2] numpy.polynomial.legendre.leggauss
[3] wiki:Gauss–Kronrod quadrature formula
[4] 高斯——勒让德公式–中南大学
[5] 各种类型的高斯积分
[6] 高斯型求积公式
[7] Calculation of Gauss Quadrature Rules
[8] Laurie D P. Calculation of Gauss-Kronrod quadrature rules[J]. Mathematics of Computation, 1997, 66(219):1133-1145
[9] Calvetti D, Golub G H, Gragg W B, et al. Computation of Gauss-Kronrod Quadrature Rules[J]. Mathematics of Computation, 2000, 69(231):1035-1052
【数值计算之二】数值积分之牛顿——科斯特公式：梯形、辛普森、辛普森3/8和布尔 高斯积分公式：勒让德、切比雪夫、拉盖尔和埃尔米特相关推荐
- 数值积分之牛顿——科斯特公式：梯形、辛普森、辛普森3/8和布尔 高斯积分公式：勒让德、切比雪夫、拉盖尔和埃尔米特
https://blog.csdn.net/sinat_21591675/article/details/86242577
- 数值计算大作业：数值积分（梯形、辛普森与龙贝格方法在Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 针对数值积分的编程,我把梯形.辛普森与龙贝格方法在MATLAB中编写的程序放在文章最后了,需要的同学自取. PS:附录中的程序并 ...
- 数值计算笔记之数值积分(一)
目录 0.引言 一.数值积分的积分思想 1.中矩形公式 2.梯形公式 3.辛普森公式 二.求积公式的余项和代数精度 三.插值型求积公式 四.牛顿--柯特斯公式 (N-C公式) 五.复化求积法 1.复化 ...
- matlab用辛普森公式求积分_数值计算实验9 数值积分实验
实验9 数值积分实验 成绩 实验类型:●验证性实验 ○综合性实验 ○设计性实验 实验目的:进一步熟练掌握变步长数值积分算法,提高编程能力和解决定积分问题的实践技能. 实验内容:用龙贝格积分算法计算 ...
- 高斯积分公式matlab_数值微分与数值积分（一）
点击蓝字 关注我们 01 数值微分与数值积分 一.数值微分 (1)数值差分与差商 微积分中,任意函数f(x)在x0点的导数是通过极限定义的: 如果去掉极限定义中h趋向于0的极限过程,得到函数在x0点处 ...
- 【数学】积分（integration）的定义，黎曼和，黎曼积分，牛顿.莱布尼茨公式，微分三大中值定理
[积分的定义] 本文将详细的描述积分的定义和牛顿.莱布尼茨公式.感谢此文:如何简单地证明.理解牛顿-莱布尼兹公式? 在上面四个图中,我们可以看到,使用等分矩形的形式来计算曲线与坐标轴所围的面积,则当细 ...
- 【数学】用C语言实现函数的定积分—— 把 “定积分定义计算出的值” 和 “牛顿-莱布尼兹公式计算出的值” 两者进行误差比较
因为考研数学看到定积分的定义以及"牛顿-莱布尼兹公式" 突然心血来潮,想用C语言把它们实现出来并对比. 1.用 "定积分定义" 计算得出数值 以及 " ...
- 数值积分： 梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式
数值积分:梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式 1.问题描述 微积分方法求积有很大的局限性,当碰到被积函数很复杂时,找不到相应的原函数.积分值 在几何上可解释为由 ...
- 数值分析复化求积matlab,MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等)...
1.理解如何在计算机上使用数值方法计算定积近似值; 2.学会复合梯形.复合Simpson和龙贝格求积分公式的编程与应用. 3.探索二重积分在矩形区域的数值积分方法. 佛山科学技术学院 实 验 报 告 ...
- 辛普森复合求积公式matlab,MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等).doc...
[摘要]佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 数值积分 专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日 一.实验目的 ...
最新文章
- UI培训分享：UI设计师要掌握哪些知识点
- markdown 创建表格
- mysql file-pos_mysql-5.7 调整mysql的复制方式由master_log_file+master_log_pos 到gtid 详解
- 问题：# mount –t ntfs /dev/sdb1 /mnt/ 解决办法
- boost::hana::is_disjoint用法的测试程序
- fastjson 过滤不需要的字段或者只要某些字段
- 计算机微机原理及接口技术实训室,《微机原理与接口技术》课程实验报告.doc...
- jQuary总结11:jQuery插件封装---jQuery封装 手风琴 动画插件
- Linux网络编程——tcp并发服务器（epoll实现）
- Java基础学习总结（91）——阿里巴巴Java开发手册公开版
- leetcode Valid Palindrome
- using在sql中是什么意思_扇贝英语地道表达法——“call for”是什么意思呢？
- 介绍几个 window 下面的terminal
- 怎么新建web程序_前端程序员发展潜力最好，那该怎么学好web前端开发？
- linux tar命令将压缩包解压到指定位置,用tar命令把目标压缩包解压到指定位置
- 电镜的成像原理-透射电镜成像原理2
- Format oracle 用法,oracle sqlplus中column格式化命令之heading用法
- 多元数量值函数积分学
- java 读取文本_Java如何读取txt文件的内容？
- m4a录音文件损坏修复_m4a录音文件怎么打开 - 卡饭网