数值积分——牛顿-科特斯公式

最新推荐文章于 2024-11-15 15:40:28 发布

努力的小流

最新推荐文章于 2024-11-15 15:40:28 发布

阅读量775

点赞数 4

分类专栏：积分与微分文章标签：数据结构

本文链接：https://blog.csdn.net/weixin_43816345/article/details/139008184

版权

积分与微分专栏收录该内容

4 篇文章

订阅专栏

本文介绍了数值积分中的牛顿-科特斯公式，包括梯形法则、复合梯形法则、辛普森1/3法则、非等距积分，并提供了在Matlab中的实现示例。此外，还讨论了Matlab的一重积分和多重积分函数的使用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

导论

关于积分的基础知识，本文不再赘述。牛顿-科特斯公式的基本思想是：使用容易积分的多项式代替复杂函数或列表数据，通用公式如下：

$I = \int_a^bf(x)dx = \int_a^bf_n(x)dx$

其中：

$f_n(x) = a_0+a_1x+...+a_{n-1}x^{n-1}+a_nx^n$

根据选取多项式的阶数，可分为梯形法则（一阶多项式）；辛普森法则（辛普森1/3法则（二次多项式）、辛普森3/8法则（三次多项式））;...

在整个区间应用一次牛顿-斯科特公式精度较低，常用的方法是将a~b积分区间分成很多小区间，在每个小区间上应用牛顿-斯科特公式，称为复合斯科特公式。

复合梯形法则

在一个区间中应用梯形法则公式如下：

$I = \int_a^b[f(a)+\frac{f(b)-f(a)}{b-a}(x-a)] = (b-a)\frac{f(a)+f(b)}{2}$

a~b之间划分n个等宽小区间，则复合梯形法则公式为：

$I = h\frac{f(x_0)+f(x_1)}{2}+h\frac{f(x_1)+f(x_2)}{2}+...+h\frac{f(x_{n-1})+f(x_n)}{2}$

$= \frac{h}{2}[f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+f(x_n)]$

下面给出一个梯形法则的Matalb实现：

function I = trap(func, a, b, n, varargin)
% 牛顿-科特斯公式求解积分（梯形法则）
if nargin<4||isempty(n),n = 100;end
x = a; h = (b-a)/n;
s = func(a, varargin{:});
for i = 1:n-1
    x = x+h;
    s = s + 2*func(x, varargin{:});
end
s = s + func(b, varargin{:});
I = (b-a) * s/(2*n);
end

复合辛普森1/3法则

当f(x0)和f(x2)之间另外添加一个中点，那么这三个点可以用抛物线相连：

$I = \int_{x_0}^{x_2}[\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\frac{(x-x_0)(x-x_1)}{(x_2-x_1)(x_2-x_0)}f(x_2)]dx$

$I = (b-a)\frac{f(x_0)+4f(x_1)+f(x_2)}{6}$

则将a~b之间划分多个等距的小区间，总的积分表示为：

$I = (b-a)\frac{f(x_0)+4\sum_{i=1,3,5}^{n-1}f(x_i)+2\sum_{j=2,4,6}^{n-2}f(x_j)+f(x_n)}{3n}$

注意：该方法中小区间的个数必须为偶数

下面给出一个辛普森1/3法则的Matlab实现：

function I = simpson(func, a, b, n, varargin)
% 辛普森1/3法则计算积分(小区间个数必须为偶数个，即n必须为偶数)
x = a; h = (b-a) / n;
s = func(a, varargin{:});
for i = 1:n-1
    x = x+h;
    if mod(i, 2) == 0
        s = s + 4*func(x, varargin{:});
    else
        s = s + 2*func(x, varargin{:});
    end
end
s = s + func(b, varargin{:});
I = (b-a) * s / (3*n);
end

辛普森3/8法则

对四个点进行三次拉格朗日插值然后积分得到：

$I = (b-a)\frac{f(x_0)+3f(x_1)+3f(x_2)+f(x_3)}{8}$

一般将辛普森1/3法则和辛普森3/8法则结合使用以应对奇数个区间的情况。

非等距积分

可使用复合梯形法则/复合辛普森法则，但无法得出最后的公式，而是通过循环在每一个小区间进行积分。

Matlab中的积分函数

一重积分

z = trapz(x, y);

其中两个向量x和y分别存储自变量和因变量，z为积分结果。

z = cumtrapz(x, y);

此函数用于计算累积积分，z为向量，其元素z(k)是从x(1)到x(k)的积分值。

trapz与cumtrapz两个函数针对离散数据点积分，而integral函数针对函数表达式进行积分。

z = integral(fun, xmin, xmax);

多重积分

q = integral2(fun, xmin, xmax, ymin, ymax);
q = integral3(fun, xmin, xmax, ymin, ymax, zmin, zmax);

q为函数fun在xmin到xmax和ymin到ymax区域上的二重积分。