MENU

Numerical Analysis: Numerical Integration

• October 7, 2019 • Read: 3946 • Knowledge,Numerical Analysis

Newton-Cotes积分

对于定积分的计算,解析解的存在特别受制于被积函数的形式,往往原函数难以找到,使得积分计算无法进行。积分的意义是对所在区域进行无穷小划分并求和,因此,从数学上来说,即使是找不到原函数,但是积分表示的物理意义依然存在。对于积分的数值方法,Newton-Cotes法是一种最容易理解的方法,这一方法实际上是一系列的该类型数值方法的总称,因为这些方法的数值积分实现的思路是一样的。

数值方法对于经典微积分的实现,一个解决的思路在于寻求近似解以替换精确解,其实在定积分的定义过程中,黎曼的定义就已经给出了这种寻求近似的思路,只不过他通过无穷小求和的形式给出了定积分的精确定义。而Newton-Cotes法就是依据这一种切割成小块然后求和的思想来实现数值解答的,其中最为关键的就是利用多项式插值来近似表示小区段内原始的函数。我们知道,这里采用的拉格朗日插值多项式,可以依据所给的插值点个数,来获得n-1次的多项式;对于2个点的区段来说,插值时只能得到一个一次多项式,这是线性插值,不管原函数在该区间上的形状是什么,都用一个一次多项式来替代;而当区间内取3个点时,拉格朗日多项式为2次,,当然也可以利用4点进行插值近似。通常,当采用两点插值的一次多项式来替换原函数时,被称为梯形法(trapz);采用三点拉格朗日插值时,被称为Simpson法;4点时称为Simpson 3/8法。

Newton-Cotes法的实现思路示意图。

trapz

trapz

所以,对于Trapz法来说,在$[x_0, x1]$这个小区间上,采用两点拉格朗日插值公式得到线性函数$f_1(x)$,那么近似解就可以表示为一次函数在这个区间上的积分。

$$ \int_{x_0} ^ {x_1} (\frac{x - x_1}{x_0 - x_1} \cdot y_0 + \frac{x - x_0}{x_1 - x_0} \cdot y_1 )\mathrm{d}x $$

因为这个积分是可以直接计算的,设区间长度(步长)为$h$,结果为:

$$ S = \frac{h}{2} \cdot y_0 + \frac{h}{2} y_1=\frac{y_0 + y_1}{2}h $$

对于这个结果,其正好时梯形的面积公式,因此当为线性插值时,该方法亦称之为梯形法

同理,采用二阶多项式时,依据拉格朗日插值多项式的要求,此时需要3个点,设为$x_0, x_1, x_2$,且取整个区段的长度为$2h$,且连个子区间长度相等,则有:

$$ \int_{x_0} ^{x_2} \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} \cdot y_0 + \frac{(x - x_0)(x - x_2)}{(x1 - x_0)(x_1 - x_2)}\ cdot y_1 + \frac{(x - x_0)(x - x_1)}{(x_2 - x_0)(x_2 - x_1)} \cdot y_2\\ $$

结果为:

$$ \frac{h}{3}(y_0 + 4y_1 + y_2) $$

这只是对一小段积分区间的计算,当每一小段的积分近似结果都得到确定后,然后进行求和,由此完成积分的数值计算。对于梯形法来说,分析的对象是由两个点组成的一个区间,当对每一段求和后,可以知道的是,整个积分区域中的第一个子区间上的$y_0$,以及最后一个区间上$y_1$均只出现一次,对于其他的点对应的$y$值来说,它们即是上一个梯形的下底,也是下一个梯形的上底,所以求和后会出现两次,这样真正在使用时的梯形法表示为:

$$ \int f(x)\mathrm{d}x = \frac{h}{2} \cdot (y_{0} + y_{n} + 2 \cdot \sum_{i = 1} ^ {n - 1} y_{i}) $$

这里,离散函数值向量元素个数为$n+1$个,index从0开始。

对于Simpson法来说,特别的一点是在研究“单个”对象时,由于其实二阶的插值,需要三个点,使得区间被分为2段,因此,就整体求解来说,该方法的子区间数目必须为偶数才能获得计算。设离散点数为$2n + 1$,则区间总的数量为$2n$。

$$ \int f(x) \mathrm{d}x = \frac{h}{3} \cdot (y_{0} + y_{2_n} + 4 \sum_{i=1} ^{n} y_{2 i -1} + 2 \cdot \sum_{i= 1} ^ {m - 1} y_{2i}) $$

这里的求和,注意区分下标为奇数以及偶数的区别。每一段中,只有中间点的函数值系数是4,并且它与前后的区段不接触,因此,单纯地保持求和就可以;除去首尾区间,区段边界都是连续接触的,所以求和会出现两次。

编程实现时,函数向量$y$的获得在外部单独计算,并且这里的计算都是等步长的。

import numpy as np

# 梯形法
def trapz_method(y, x):
    h = x[1] - x[0]
    s = h / 2 * (y[0] + y[-1] + 2 * y[1:-1].sum())
    return s

# 以x**2在[0, 1]上的积分为例
x = np.linspace(0, 1, 201)
y = x**2
res = trapz_method(y, x)
'''Simpson法中必须保证离散的点个数为奇数,才能确保
区间段数为偶数。
'''

def simpson_method(y, x):
    h = x[1] - x[0]
    s = h / 3 * (y[0] + 
                 y[-1] + 
                 4 * y[1:-1:2].sum() +
                 2 * y[2:-1:2].sum())
    return s
Archives Tip
QR Code for this page
Tipping QR Code