一、数学原理
二、实验内容
本实验由1个小实验组成,在进行积分计算时,若需要精度,可取ε=1e-8.积分的计算结果保留六位有效数字。
分别取5个节点和8个节点,用牛顿-柯特斯公式计算积分
代码实现:
function y=NewtonCotes(fun,a,b,n)
if nargin==1
[mm,nn]=size(fun);
if mm>=8
error('为了保证NewtonCotes积分的稳定性,最多只能有9个等距节点!')
elseif nn~=2
error('fun构成应为:第一列为x,第二列为y,并且个数为小于10的等距节点!')
end
xk=fun(1,:);
fk=fun(2,:);
a=min(xk);
b=max(xk);
n=mm-1;
elseif nargin==4
% 计算积分节点xk和节点函数值fx
xk=linspace(a,b,n+1);
if isa(fun,'function_handle')
for i=1:n+1
fx(i)=fun(xk(i));
end
else
error('fun积分函数的句柄')
end
else
error('输入参数错误,请参考函数帮助!')
end
% 计算柯特斯系数
format rat %以分数形式输出
Ck=cotescoeff(n);
% 计算求积系数
Ak=(b-a)*Ck;
% 求和算积分
y=Ak*fx';
format long
function Ck=cotescoeff(n)
%{
% 由于柯特斯系数最多7阶,可以直接使用
switch n
case 1
Ck=[1,1]/2;
case 2
Ck=[1,4,1]/6;
case 3
Ck=[1,3,3,1]/8;
case 4
Ck=[7,32,12,32,7]/90;
case 5
Ck=[19,75,50,50,75,19]/288;
case 6
Ck=[41,216,27,272,27,216,41]/840;
case 7
Ck=[751,3577,1323,2989,2989,1323,3577,751]/17280;
otherwise
fprintf('n=%d',n), error('错误的柯特斯阶数')
end
%}
% 但是为了体现公式2.2,也可以使用程序计算n阶柯特斯系数
for i=1:n+1
k=i-1;
Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n);
end
function f=intfun(t,n,k)
% 柯特斯系数中的积分表达式
f=1;
for i=[0:k-1,k+1:n]
f=f.*(t-i);
end
例题:
Example
format compact %(有3处需要根据注释填写代码)
clc,clear
ep=1e-8;
fun=@ (x)sin(x.*x);%(1)填写函数形式,必须可以接受矢量输入
y1= NewtonCotes(fun,0,1,4); %(2)5节点
y2= NewtonCotes(fun,0,1,7); %(3)8节点
%若无结果,请仔细阅读指导书的实验分析1,再修改代码使其有结果
y3=integral(fun,0,1); %matlab提供的求积分函数,可认为是相对精确的
fprintf('牛顿-柯特斯求积结果,5节点:%.6f\n',y1);
fprintf('牛顿-柯特斯求积结果,8节点:%.6f\n',y2);
fprintf(' 较精确值:%.6f\n',y3);
结果与分析: