数值积分——(Newton-Cotes)牛顿-科特斯求积公式 | 北太天元

一、Newton-Cotes求积公式

对 $[a,b]$ 的 $n$ 等分点 $x_k=a+kh,h=\frac{b-a}{n},k=0,1,2,\cdots ,n$

$n$ 阶Newton-Cotes 公式
$${\int }_a^{b}f(x)\mathrm{d}x\approx (b-a)\sum _{k=0}^n{c}_k^{(n)}f(x_k)$$

Cotes 系数，令$x=a+th$
$$C_k^{(n)}=\frac{1}{n}{\int }_0^n\left(\prod _{\genfrac{}{}{0ex}{}{i=0}{i\ne k}}^n\frac{t-i}{k-i}\right)\mathrm{d}t=\frac{(-1{)}^{n-k}}{nk!(n-k)!}{\int }_0^n\prod _{\stackrel{i=0}{i\ne k}}^n(t-i)\mathrm{d}t$$

Cotes系数表

n	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$	$C_k^{(n)}$
1	$\frac{1}{2}$	$\frac{1}{2}$
2	$\frac{1}{6}$	$\frac{4}{6}$	$\frac{1}{6}$
3	$\frac{1}{8}$	$\frac{3}{8}$	$\frac{3}{8}$	$\frac{1}{8}$
4	$\frac{7}{90}$	$\frac{32}{90}$	$\frac{2}{90}$	$\frac{32}{90}$	$\frac{7}{90}$
5	$\frac{19}{228}$	$\frac{75}{228}$	$\frac{50}{228}$	$\frac{50}{228}$	$\frac{75}{228}$	$\frac{19}{228}$
6	$\frac{41}{840}$	$\frac{216}{840}$	$\frac{27}{840}$	$\frac{272}{840}$	$\frac{27}{840}$	$\frac{216}{840}$	$\frac{41}{840}$
7	$\frac{751}{17280}$	$\frac{3577}{17280}$	$\frac{1323}{17280}$	$\frac{2989}{17280}$	$\frac{2989}{17280}$	$\frac{1323}{17280}$	$\frac{3577}{17280}$	$\frac{751}{17280}$
8	$\frac{989}{28350}$	$\frac{5888}{28350}$	$-\frac{928}{28350}$	$\frac{10496}{28350}$	$-\frac{4540}{28350}$	$\frac{10496}{28350}$	$-\frac{928}{28350}$	$\frac{5888}{28350}$	$\frac{989}{28350}$
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$

n=1 时，代入Cotes系数得到梯形公式
$${\int }_a^{b}f(x)\mathrm{d}x\approx \frac{b-a}{2}(f(a)+f(b))$$
n=2 时,代入Cotes系数得到 Simpson公式
$${\int }_a^{b}f(x)\mathrm{d}x\approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right).$$
n=4 时,代入Cotes系数得到 四阶Newton-Cotes公式
$${\int }_a^{b}f(x)\mathrm{d}x\approx \frac{b-a}{90}(7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)).$$

二、算法

$♡$ Newton-Cotes求积公式：[NC] = ncotes_integral(a,b,n,f)

1. 输入

   • $[a,b]$
   • $n$：将$[a,b]$ $n$ 等分
   • $f$：已经定义好的函数，支持向量运算

2. 实现步骤

   • 通过$[a,b]$ $n$ 等分获得 $n+1$ 个横坐标构成的$x_0$向量
   • $y_0=f(x_0)$
   • $n$ 确定 Newton-Cotes求积公式的阶数，和所要用到的 Cotes 系数
   • 代入
     $${\int }_a^{b}f(x)\mathrm{d}x\approx (b-a)\sum _{k=0}^n{c}_k^{(n)}f(x_k)=NC$$

3. 输出

   • $NC$：通过Newton-Cotes公式得到的积分近似值

三、北太天元源程序

function [NC] = ncotes_integral(a,b,n,f)
% Newton-Cotes 求积公式
% [a,b] 的 n 等分
% f：提前定义好的函数,要求支持向量运算
% n 不超过 8
%linspace可以把[a,b]等分成n个点，n-1个区间
%
%   Version:            1.0
%   last modified:      07/06/2023
x0 = linspace(a,b,n+1); % 故此处是n+1
y0 = f(x0);
Cotes = cell(1,8); % 创建一个空的元胞数组
Cotes{1} = [1/2 1/2];
Cotes{2} = [1/6 4/6 1/6];
Cotes{3} = [1/8 3/8 3/8 1/8];
Cotes{4} = [7/90 32/90 2/90 32/90 7/90];
Cotes{5} = [19/228 75/228 50/228 50/228 75/228 19/228];
Cotes{6} = [41/840 216/840 27/840 272/840 27/840 216/840 41/840];
Cotes{7} = [751/17280 3577/17280 1323/17280 2989/17280 2989/17280 1323/17280 3588/17280 751/17280];
Cotes{8} = [989/28350 5888/28350 -928/28350 10496/28350 -4540/28350 10496/28350 -928/28350 5888/28350 989/28350];

sum_yc = sum(Cotes{n} .* y0);
NC = (b-a) * sum_yc;
end

保存为ncotes_integral.m文件

四、数值算例

用Newton-Cotes公式计算
$${\int }_{-4}^4\frac{\mathrm{d}x}{1+x^2}=2\arctan (4)\approx 2.6516$$

clc,clear all,format long;
f1 = @(x)1./(1+x.^2);
a = -4; b = 4;

real = 2.6516; % 真实值取4位小数的值
Nc = zeros(1,8);
delta = zeros(1,8);
for n =1:1:8
    Nc(n) = ncotes_integral(a,b,n,f1);
    delta(n) = abs(Nc(n)-real);
end
n =1:8;
figure(1);
plot(n,Nc,'-*b');
hold on
    plot(n,real*ones(1,8),'-o')
legend('N-C','real')
title('不同n下N-C与真实值的区别')

figure(2);
plot(n,delta,'-*r');

运行后得


可以看到 n=8 时， 误差不降反增，这主要是因为Cotes系数中出现了负数，使其稳定性下降

实际应用时，只会使用 $n\le 7$ 的求积公式。