Python数值计算（34）——newton-cotes积分公式

C囧囧

已于 2025-01-26 17:56:42 修改

阅读量417

点赞数 8

分类专栏： Python 文章标签： python 开发语言数值积分

于 2024-10-29 22:09:49 首次发布

本文链接：https://blog.youkuaiyun.com/cdinten/article/details/143327177

版权

Python 专栏收录该内容

45 篇文章

订阅专栏

1. 背景知识

前面可以看到，通过使用2个到4个不同的等距节点，采用Lagrange插值公式构造多项式曲线拟合原待求积分的函数，实现数值积分，从两点的梯形公式，到三点的1/3 Simpson公式，以及四点的3/8 Simpson公式，其本质都是一样的，按这种思路，我们可以一直做下去，例如五点、六点等的积分公式，那么，这其中的所蕴含的数学思想，就是Newton-Cotes公式。

首先假设在 $[a,b]$ 上有 $(n+1)$ 个等距点，即：

$x_k=a+kh,h=\frac{b-a}{n},0\leq k\leqslant n$

由前面介绍的曲线拟合中的插值算法，可知这些点及其对应函数值，通常可以构成次数不超过n次的多项式，令 $x=a+th$ 进行换元，其中t为变量，则可得：

$\int_{a}^{b}l_k(x)dx=\int_{a}^{b}(\prod_{j=0,j\neq k}^{n}\frac{x-x_j}{x_k-x_j})dx\\ =\int_{0}^{n}(\prod_{j=0,j\neq k}^{n}\frac{t-j}{k-j})hdt\\ =\frac{(-1)^{n-k}h}{k!(n-k)!}\int_{0}^{n} \prod_{j=0,j\neq k}^{n}(t-j)dt\\ =(b-a)\frac{(-1)^{n-k}}{k!(n-k)!n}\int_{0}^{n} \prod_{j=0,j\neq k}^{n}(t-j)dt\\ =(b-a)*c_{k}^{(n)}$

其中 $c_{k}^{(n)}$ 就是Newton-Cotes求积系数，改值只与n和k有关，而与被积函数，积分区间均无关。

2. 系数计算

使用前面介绍的算法，计算不同阶次系数的代码如下：

from fractions import Fraction
from numpy.polynomial.polynomial import Polynomial

def factorial(n):
    if n==0:
        return 1
    else:
        p=1
        for i in range(1,n+1):
            p*=i
        return p

def nc_coef(n, k):
    if k<n/2:
        return nc_coef(n,n-k)
    p=Polynomial([1])
    for j in range(0,n+1):
        if j==k:
            continue
        p=p*Polynomial([-j,1])
    fn1=factorial(n+1)
    pInt=p.integ()*fn1
    ff=round(pInt(n)-pInt(0))
    f=Fraction(int(pInt(n)-pInt(0)),fn1)
    f=f/factorial(k)/factorial(n-k)/n
    if (n-k)%2==1:
        f=-f
    return f


for n in range(1,10):
    for k in range(0,n+1):
        print(nc_coef(n,k),end=' ')
    print()

输出为：

1/2 1/2 
1/6 2/3 1/6
1/8 3/8 3/8 1/8
7/90 16/45 2/15 16/45 7/90
19/288 25/96 25/144 25/144 25/96 19/288 
41/840 9/35 9/280 34/105 9/280 9/35 41/840
751/17280 3577/17280 49/640 2989/17280 2989/17280 49/640 3577/17280 751/17280 
989/28350 2944/14175 -464/14175 5248/14175 -454/2835 5248/14175 -464/14175 2944/14175 989/28350
2857/89600 15741/89600 27/2240 1209/5600 2889/44800 2889/44800 1209/5600 27/2240 15741/89600 2857/89600

可见，

$n=1$ 时，是两点确定的一阶多项式，这个时候就是梯形公式；

$n=2$ 时，是三点确定的二阶多项式，这个时候就是Simpson 1/3公式；

$n=3$ 时，是四点确定的三阶多项式，这个时候就是Simpson 3/8公式；

$n=4$ 时，虽然前面没有介绍，但是这个被称为Cotes公式。

3. 演示

以下展示了n在不同取值时，对humps函数在[0,1]上进行数值积分的效果：

实现源代码如下：

from numpy.polynomial import Polynomial
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from humps import humps,AnalyticalIntOfHumps

def genpoly(f,a,b,n):
    x=np.linspace(a,b,n)
    y=f(x)
    p=Polynomial.fit(x,y,n-1)
    return p

fig,ax=plt.subplots()
plt.subplots_adjust(bottom=0.25)
ax_n=plt.axes([0.15,0.1,0.7,0.03])
s_n=Slider(ax_n,'n',2,10,valinit=3,valstep=1)
a=0
b=1
def update(val):
    n=int(s_n.val)
    p=genpoly(humps,a,b,n)
    x=np.linspace(a,b,1000)
    y=p(x)
    ax.clear()
    ax.plot(x,y,'g--')
    ax.fill_between(x,0,y,color='lime',alpha=0.62)

    ax.plot(x,humps(x),'r-')

    t=np.linspace(a,b,n)
    ax.plot(t,humps(t),'b.')

    anaInt=AnalyticalIntOfHumps(a,b)
    polyInt=p.integ()(b)-p.integ()(a)
    ax.text(1,70,f'anaInt={anaInt:.6f}\npolyInt={polyInt:.6f}\nrelErr={abs(polyInt-anaInt)/anaInt*100:.6f}%',va='bottom',ha='right')
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    fig.canvas.draw_idle()

s_n.on_changed(update)
update(0)
plt.show()