牛顿-科特斯(Newton-Cotes)数值积分学习笔记

微积分是我们数学领域一个强大的工具，科学研究领域乃至生活中的方方面面都渗透着微积分的思想。
学过高等数学的我们都知道，对于$I={\int }_a^{b}f(x)$这样一个定积分，运用牛顿-莱布尼兹公式（Newton-Leibniz formula)$${\int }_a^{b}f(x)=F(b)-F(a)$$可以很方便快捷地求解一些较为简单的连续性问题。但现实中，大多数问题都是离散地、非连续地，我们没法通过代数手段找到被积函数所对应的原函数。因此，$数值积分$成了我们求解积分问题的一个有力手段。
积分中值定理大家想必很熟悉了：若函数$f(x)$在闭区间$[a,b]$上连续，则在积分区间上至少存在一个点$\xi$，使下式成立
$${\int }_a^{b}f(x)=(b-a)f(\xi )$$

这个定理的几何意义便是：$存在一个矩形，其面积与连续曲线与坐标轴所围成的封闭图形的面积相等$。
也就是说，当我们能找到这一点$\xi$时及其对应的$f(\xi )$时，被积函数在一段区间上的积分值就转化成了某一点的函数值与一个定值常数的乘积。这样一来，原来的积分仿佛被一种数值形式所取代，数值积分便由此有了雏形。这便是所谓的离散化思想。 相应地，可以简单定义这种普遍存在的数值形式：$${\int }_a^{b}f(x)=A_{i}f(x_i)$$
这里，$f(x_i)$表示某点函数值,$A_i$代表该函数值对应权重。

然而，$\xi$是多少，$f(\xi )$又是多少，我们现在凭借以上并不能确定，于是牛顿-科特斯(Newton-Cotes)公式将逐步求解出我们的这些疑惑。
首先思考的问题便是怎么求解出这些所谓的权重，观察上述式子，注意到左边为积分项，右边为数值项，考虑是否也能将左半部分的积分项也化成数值项的形式，我们将用到拉格朗日插值（Lagrange Interpolation Polynomial）。

拉格朗日插值

对于某个未知的函数，已知该函数上一系列未知的散点$\left(x_0,y_0\right),\left(x_1,x_1\right),\cdots ,\left(x_n,y_n\right)$
定义基函数：$$l_i(x)=\frac{(x-x_0)\cdots (x-x_{i-1})(x-x_{i+1})\cdots (x-x_n)}{(x_i-x_0)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}$$ $$=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j},j=0,1,2,\cdots ,n$$
该函数所满足的就是在散点$x=x_i$处，$l_i(x)=1$,而其他点$x_j$处均为$0$.
于是，拉格朗日插值函数：
$L_n(x)={\sum }_{i=0}^n{y}_i{l}_i(x)$
从而巧妙地用函数的形式刻画了这种$曲线经过固定点$现象。
这里，参考MATLAB之拉格朗日插值

function y=Lagrange1(x0,y0)
syms x

for i=1:length(x0)
  x0_=x0;
    x0_(i)=[];
    a{i}=-x0_+x0(i);
    b=-x0_+x;
    l(i)=prod(b)/prod(a{i}); x0_=x0;
end
for i=1:length(x0)
    L(i)=y0(i)*l(i);
end
Ln=sym2poly(expand(sum(L)));
xx=x0(1):0.1:x0(end);
y=polyval(Ln,xx);
Y1=interp1(x0,y0,xx,'spline');
Y2=interp1(x0,y0,xx,'nearest');
plot(x0,y0,'r*','LineWidth',2);hold on
plot(xx,Y1,'LineWidth',2);
plot(xx,y,'LineWidth',2);
plot(xx,Y2,'LineWidth',2);hold off
xlabel('x');ylabel('y');
legend('data','spline','Lagrange Interpolation Polynomial','nearest')

这里我们进行了对同一组数据点拉格朗日插值和最近项插值以及样条线插值的绘制，可以发现拉格朗日插值与样条曲线已大致合为一体。

有了以上的知识储备，我们可以对$f(x)$进行插值运算，化其为
离散多项式，即：
$$f(x)=L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)$$
于是，$${\int }_a^{b}f(x)dx={\int }_a^b\sum _{i=0}^n{y}_i{l}_i(x)dx={\int }_a^b\sum _{i=0}^n{f}_i(x)l_i(x)dx=\sum _{i=0}^n{f}_i(x){\int }_a^b{l}_i(x)dx$$
记$${\int }_a^b{l}_i(x)dx=A_i={\int }_a^b\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}dx\text{\hspace{1em}(1)}$$
这样求积公式便由最初的$${\int }_a^{b}f(x)=A_{i}f(x_i)$$转化成了更具有一般规律的
$${\int }_a^{b}f(x)=\sum _{i=0}^n{A}_{i}f(x_i)$$
观察上述$(1)$式，任意未知函数的积分形式将演化为幂函数的积分形式，无论从那种角度来说，幂函数的积分将大大简化积分繁琐性。
根据之前的拉格朗日插值函数，我们可在此基础上求解牛顿-科特斯(Newton-Cotes)数值积分。

function y=Newton_Cotes(x0,y0,qa,qb)
syms x
for i=1:length(x0)
  x0_=x0;
    x0_(i)=[];
    a{i}=-x0_+x0(i);
    b=-x0_+x;
    l(i)=prod(b)/prod(a{i});
    A(i)=int(l(i),qa,qb);
    x0_=x0;
end
%这里用于定义步长，不同步长将有不同求解精度
% 步长=(x(end)-x(1))/1-------一阶---------梯形公式
%步长=(x(end)-x(1))/2-------二阶----------辛普森求解公式
%步长=(x(end)-x(1))/4-------四阶----------科斯特公式

for i=1:length(x0)
    I(i)=y0(i)*A(i);
end
y=eval(sum(I));

我们用以上定义的$functiony=NewtonCotes(x0,y0,qa,qb)$分别测试了一次函数、二次函数、正弦均得到问题的精确值。
二次抛物线：$$y=x^2$$
列些曲线上的任意点$x0=[1,2,3]，y0=[1,4,9]$
MATLAB求解$${\int }_1^4{x}^2=\sum _{i=0}^n{A}_{i}f(x_i)=NewtonCotes(x0,y0,1,4)$$
计算结果为21为该问题的精确解答。

考虑到$(1)$幂函数积分项是否也能化为幂函数之和，于是将$[a,b]$区间划分成$n$等份，每份步长为$h=\frac{b-a}{2}$,数值积分点取$x_i=a+kh,k=0,1,2,\cdots ,n$
在$${\int }_a^b{l}_i(x)dx=A_i={\int }_a^b\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}dx\text{\hspace{1em}(1)}$$中令$x=a+th$,于是：
$$A_i=\frac{b-a}{n}\frac{(-1{)}^{n-i}}{i!(n-i)!}{\int }_0^n\prod _{j=0,j\ne i}^n(t-j)dt,i=0,1,2,\cdots ,n$$
记$$C_i=\frac{(-1{)}^{n-i}}{ni!(n-i)!}{\int }_0^n\prod _{j=0,j\ne i}^n(t-j)dt,i=0,1,2,\cdots ,n$$为科斯特系数，显然科斯特系数与被积函数与积分区间均无关。

$n$取1，2，4，该数值求积公式将分别转化成$梯形公式、辛普森公式、科斯特公式$。

得到的科斯特系数表：