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

清水折木

已于 2024-07-03 23:04:37 修改

阅读量2.6k

点赞数 28

分类专栏：数值计算方法 | 北太天元 or Matlab 文章标签：线性代数数值计算数值积分

于 2024-04-25 11:01:15 首次发布

本文链接：https://blog.csdn.net/Math_Boy_W/article/details/138181698

版权

数值计算方法 | 北太天元 or Matlab 专栏收录该内容

17 篇文章 9 订阅

订阅专栏

本文介绍了Newton-Cotes求积公式，包括其基本原理、Cotes系数表以及Matlab中的算法实现。通过实例展示了如何使用这些公式计算积分并分析其在不同阶数下的精度。

摘要由CSDN通过智能技术生成

文章对应视频讲解：(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^bf(x)\mathrm{d}x\approx(b-a)\sum_{k=0}^nc_k^{(n)}f(x_k)$

Cotes 系数，令 $x = a + t h$
$C_{k}^{(n)}={\frac{1}{n}}\int_{0}^{n}\left(\prod_{i=0\atop i\neq k}^{n}{\frac{t-i}{k-i}}\right)\mathrm{d}t=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_{0}^{n}\prod_{\overset{i=0}{i\neq 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	$\dfrac{1}{2}$	$\dfrac{1}{2}$
2	$\dfrac{1}{6}$	$\dfrac{4}{6}$	$\dfrac{1}{6}$
3	$\dfrac{1}{8}$	$\dfrac{3}{8}$	$\dfrac{3}{8}$	$\dfrac{1}{8}$
4	$\dfrac{7}{90}$	$\dfrac{32}{90}$	$\dfrac{12}{90}$	$\dfrac{32}{90}$	$\dfrac{7}{90}$
5	$\dfrac{19}{288}$	$\dfrac{75}{288}$	$\dfrac{50}{288}$	$\dfrac{50}{288}$	$\dfrac{75}{288}$	$\dfrac{19}{288}$
6	$\dfrac{41}{840}$	$\dfrac{216}{840}$	$\dfrac{27}{840}$	$\dfrac{272}{840}$	$\dfrac{27}{840}$	$\dfrac{216}{840}$	$\dfrac{41}{840}$
7	$\dfrac{751}{17280}$	$\dfrac{3577}{17280}$	$\dfrac{1323}{17280}$	$\dfrac{2989}{17280}$	$\dfrac{2989}{17280}$	$\dfrac{1323}{17280}$	$\dfrac{3577}{17280}$	$\dfrac{751}{17280}$
8	$\dfrac{989}{28350}$	$\dfrac{5888}{28350}$	$-\dfrac{928}{28350}$	$\dfrac{10496}{28350}$	$-\dfrac{4540}{28350}$	$\dfrac{10496}{28350}$	$-\dfrac{928}{28350}$	$\dfrac{5888}{28350}$	$\dfrac{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^bf(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^bf(x)\mathrm{d}x\approx\frac{b-a}{90}(7f(x_0)+32f(x_1)+12f(x_2)+32f(x_3)+7f(x_4)).$

二、算法

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

输入
- $[a, b]$
- $n$ ：将 $[a, b]$ $n$ 等分
- $f$ ：已经定义好的函数，支持向量运算
实现步骤
- 通过 $[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$
输出
- $NC$ ：通过Newton-Cotes公式得到的积分近似值

三、北太天元 or matlab 实现

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 12/90 32/90 7/90];
Cotes{5} = [19/288 75/288 50/288 50/288 75/288 19/288];
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)\approx2.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');