复化科特斯公式matlab_【原创】牛顿-柯特斯数值积分公式及其MATLAB的实现

一、数值积分基本公式

数值求积基本通用公式如下

Eqn1.gif (1.63 KB, 下载次数: 40)

2009-11-20 23:23 上传

xk：求积节点

Ak：求积系数，与f(x)无关

数值积分要做的就是确定上式中的节点xk和系数Ak。可以证明当求积系数Ak全为正时，上述数值积分计算过程是稳定。

二、插值型数值积分公式

对f(x)给定的n+1个节点进行Lagrange多项式插值，故

Eqn2.gif (2.95 KB, 下载次数: 32)

2009-11-20 23:23 上传

即求积系数为

Eqn3.gif (3.29 KB, 下载次数: 32)

2009-11-20 23:23 上传

三、牛顿-柯特斯数值积分公式

当求积节点在[a,b]等间距分布时，插值型积分公式(先使用Lagrange对节点进行多项式插值，再计算求积系数，最后求积分值)称为Newton-Cotes积分公式。

由于Newton-Cotes积分是通过Lagrange多项式插值变化而来的，我们都知道高次多项式插值会出现Runge振荡现象，因此会导致高阶Newton-Cotes公式不稳定。

Newton-Cotes积分公式的求积系数为

Eqn4.gif (3.38 KB, 下载次数: 31)

2009-11-20 23:28 上传

其中C(k,n)称为柯特斯系数。

(1)当n=1时，Newton-Cotes公式即为梯形公式

Eqn5.gif (1.68 KB, 下载次数: 34)

2009-11-20 23:28 上传

容易证明上式具有一次代数精度(对于Newton-Cotes积分公式，n为奇数时有n次迭代精度，n为偶数时具有n+1次精度，精度越高积分越精确，同时计算量也越大)

(2)当n=2时，Newton-Cotes公式即为辛普森(Simpson)公式或者抛物线公式

Eqn6.gif (2.04 KB, 下载次数: 30)

2009-11-20 23:28 上传

上式具有3次迭代精度

(3)当n=4时，Newton-Cotes公式称为科特斯(Cotes)公式

Eqn7.gif (2.68 KB, 下载次数: 30)

2009-11-20 23:28 上传

上式具有5次迭代精度。由于n=3和n=2时具有相同的迭代精度，但是n=2时计算量小，故n=3的Newton-Cotes积分公式用的很少

(4)当≥8时，通过计算可以知道，在n=8时柯特斯系数出现负值

由于数值积分稳定的条件是求积系数Ak必须为正，所以n>=8以上高阶Newton-Cotes公式，我们不能保证积分的稳定性(其根本原因是，Newton-Cotes公式是由Lagrange插值多项推导出来的，而高阶多项式会出现Rung现象)。

四、复化求解公式

n阶Newton-Cotes公式只能有n+1个积分节点，但是高阶Newton-Cotes公式由不稳定。为了提高大区间的数值积分精度，我们采用了分段积分的方法，即先将原区间划分成若干小区间，然后对每一个小区间使用Newton-Cotes积分公式，这就是复化Newton-Cotes求积公式。

(1)当n=1时，称为复化梯形公式。将[a,b]等分为n份，子区间长度为h=(b-a)/n，则复化梯形公式为

(注意：复化求解公式不需要求积子区间等间距，只是Newton-Cotes公式分段积分时自动对小区间进行等分，我们这里采用等分子区间是为了便于计算而已)

Eqn8.gif (2.18 KB, 下载次数: 28)

2009-11-20 23:28 上传

(2)当n=2时，称为复化辛普森公式。

Eqn9.gif (2.96 KB, 下载次数: 25)

2009-11-20 23:28 上传

五、Newton-Cotes数值积分公式MATLAB代码

复化Newton-Cotes数值积分公式

function y=mulNewtonCotes(fun,a,b,m,n)

% 复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和

% 参数说明

% fun，积分函数的句柄，必须能够接受矢量输入

% a，积分下限

% b，积分上限

% m，将区间[a,b]等分的子区间数量

% n，采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性

%    (1)n=1，即复化梯形公式

%    (2)n=2，即复化辛普森公式

%    (3)n=4，即复化科特斯公式

%

% Example

% fun=@(x)sin(x).*cos(x)

% mulNewtonCotes(fun,0,2,10,4)

%

% by dynamic of Matlab技术论坛

% see also http://www.matlabsky.com

% contact me matlabsky@gmail.com

% 2009-11-20 22:35:32

%

xk=linspace(a,b,m+1);

for i=1:m

s(i)=NewtonCotes(fun,xk(i),xk(i+1),n);

end

y=sum(s);复制代码牛顿-科特斯数值积分公式(全部代码参见附件)

function [y,Ck,Ak]=NewtonCotes(fun,a,b,n)

% y=NewtonCotes(fun,a,b,n)

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

%

% 参数说明：

% fun，积分表达式，这里有两种选择

%      (1)积分函数句柄，必须能够接受矢量输入，比如fun=@(x)sin(x).*cos(x)

%      (2)x,y坐标的离散点，第一列为x，第二列为y，必须等距，且节点的个数小于9，比如：fun=[1:8;sin(1:8)]'

% 如果fun的表采用第二种方式，那么只需要输入第一个参数即可，否则还要输入a,b,n三个参数

% a，积分下限

% b，积分上限

% n，牛顿-科特斯数公式的阶数，必须满足1<=n<=7，因为n>=8时不能保证公式的稳定性

%    (1)n=1，即梯形公式

%    (2)n=2，即辛普森公式

%    (3)n=4，即科特斯公式

% y，数值积分结果

% Ck，科特斯系数

% Ak，求积系数

%

% Example

% fun1=@(x)sin(x);%必须可以接受矢量输入

% fun2=[0:0.1:0.5;sin(0:0.1:0.5)];%最多8个点，必须等距

% y1=NewtonCotes(fun1,0,0.5,6)

% y2==NewtonCotes(fun2)

%

% by dynamic of Matlab技术论坛

% see also http://www.matlabsky.com

% contact me matlabsky@gmail.com

% 2009-11-20 15:06:51复制代码