牛顿.docx
《牛顿.docx》由会员分享，可在线阅读，更多相关《牛顿.docx（11页珍藏版）》请在冰豆网上搜索。
牛顿
牛顿-柯特斯数值积分公式及其MATLAB的实现
（2009-11-2023:
46:
18）
至于Newton-Cotes数值积分的具体表达、推导、说明和附件，下面给只是给出Newton-Cotes数值积分的Matlab实现代码
复化Newton-Cotes数值积分公式
functiony=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）
xk=linspace（a,b,m+1）;
fori=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=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）
%
%bydynamicofMatlab技术论坛
%seealso
%contactme matlabsky@
%2009-11-2015:
06:
51
%
ifnargin==1
[mm,nn]=size（fun）;
ifmm>=8
error（'为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！
'）
elseifnn~=2
error（'fun构成应为：
第一列为x，第二列为y，并且个数为小于10的等距节点！
'）
end
xk=fun（1,:
）;
fk=fun（2,:
）;
a=min（xk）;
b=max（xk）;
n=mm-1;
elseifnargin==4
%计算积分节点xk和节点函数值fx
xk=linspace（a,b,n+1）;
ifisa（fun,'function_handle'）
fx=fun（xk）;
else
error（'fun积分函数的句柄，且必须能够接受矢量输入！
'）
end
else
error（'输入参数错误，请参考函数帮助！
'）
end
%计算科特斯系数
Ck=cotescoeff（n）;
%计算求积系数
Ak=（b-a）*Ck;
%求和算积分
y=Ak*fx';
functionCk=cotescoeff（n）
%由于科特斯系数最多7阶，为了方便我们可以直接使用，省得每次都计算
%A1=[1,1]/2
%A2=[1,4,1]/6
%A3=[1,3,3,1]/8
%A4=[7,32,12,32,1]/90
%A5=[19,75,50,50,75,19]/288
%A6=[41,216,27,272,27,216,41]/840
%A7=[751,3577,1323,2989,2989,1323,3577,751]/17280
%当时为了体现公式，我们使用程序计算n阶科特斯系数
fori=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
functionf=intfun（t,n,k）
%科特斯系数中的积分表达式
f=1;
fori=[0:
k-1,k+1:
n]
f=f.*（t-i）;
end
求f（x）在[a,b]上的定积分用柯特斯公式
把[a,b]4等分步长h=（b-a）/4取等分点
x[i]=a+i*h（i=0,1,2,3,4）
柯特斯公式为
C=1/90*（b-a）*（7*f（x[0]）+32*f（x[1]）+12*f（x[2]）+32*f（x[3]）+7*f（x[4]））
复化求积分
先将[a,b]等分成n份步长为h=（b-a）/n,x[k]=a+k*h（k=0,1,2...n）
先用柯特斯公式在每个子段[x[k],x[k+1]]求得积分c[k]
然后求数列c[0],c[1],c[2]...c[n-1]的和作为积分
∙
∙
∙zhulinpptor
∙（zhulin）
∙等 级：
∙
#
#include
#include
typedefdouble（*pFun）（double）;
doublefun（doublex）
{
returnx*x*x*x*x*x;
}
doubleCotes（doublea,doubleb）//柯特斯公式有5次代数精度
//如果机械求积公式对x^k均能准确成立，
//但对不准确x^（k+1）则称机械求积公式具有k次代数精度
{
doublec[5];doublevalue=0.0;
intd[5]={7,32,12,32,7};
doubleh=（b-a）/4.0;
for（inti=0;i<5;i++）{c[i]=a+i*h;}
pFunf=fun;
for（i=0;i<5;i++）{value+=f（c[i]）*d[i];}
returnvalue*（b-a）/90.0;
}
doubleComplex（doublea,doubleb,doublee）//e为误差//复合柯特斯公式求积分
{
doublep1,p2;
intn=2;inti;
doubleh;
p1=Cotes（a,b）;
p2=Cotes（a,（a+b）/2）+Cotes（（a+b）/2,b）;
while（fabs（p1-p2）>e）
{
p1=p2;
n=2*n;h=（b-a）/n;p2=0;
for（i=0;i}
returnp2;
}
intmain（）
{
doublearea;
area=Cotes（0,8）;
printf（"Cotes:
area=%f\n",area）;
area=Complex（0,8,0.0001）;
printf（"Complex:
area=%f\n",area）;
area=Complex（0,8,0.00000001）;
printf（"Complex:
area=%f\n",area）;
return1;
}
C++编程用梯形求积公式求解定积分∫3lnxdx积分区间为（1，2）的值
#include"stdio.h"
#include"math.h"
#defineN100
voidmain（）
{
doubledelta=1.0/N;
doublex=1;
doubley1,y2;
doubles=0;
y1=0;
while（x<2）
{
y2=3*log（x+delta）;
s+=（y1+y2）*delta/2;
y1=y2;
x+=delta;
}
printf（"%lf\n",s）;
}
#include
#include
usingnamespacestd;
intmain（）{
doublex=1;
doubleh=1e-6;
doublesum=0;
while（x<2）{
sum+=log（x）;
x+=h;
}
sum=sum*h*3;
cout<return0;
}
MATLAB利用复合梯形公式求解积分
可以利用matlab的trapz函数命令
x=0:
0.00001:
1;%x用来储存积分点
y=（x+1）.*sin（x）;%y用来求解积分点x处的函数值
I=trapz（x,y）
I=
0.7608663730793
验证该问题的解析解
symsx
y=（x+1）*sin（x）;%被积函数表达式
II=int（y,0,1）
II=
sin
（1）-2*cos
（1）+1%II即为该被积函数的解析解
II_E=eval（II）
II_E=
0.760866373071617%II的数值解
%可以看出梯形求积公式在步长等于0.00001的情况下，数值积分的解与解析解的数值能达到小数点后11位保持一致
数值能达到小数点后11位保持一致
赞同
>>chazhijifen
x=127.5040
tixing
积分精确值Iexact=3.141592653
nIError
23.1000000000.041592653
43.1311764710.010416182
83.1389884940.002604159
163.1409416120.000651041
323.1414298930.000162760
643.1415519630.000040690
1283.1415824810.000010172
2563.1415901100.000002543
5123.1415920180.000000635
simpson
Simpson法计算积分
精确值Iexact=3.141592653;
nIError
23.1333333330.008259320
43.1415686270.000024026
83.1415925020.000000151
163.1415926510.000000002
323.141592654-0.000000001
643.141592654-0.000000001
1283.141592654-0.000000001
2563.141592654-0.000000001
5123.141592654-0.000000001
复化辛甫生求积公式（Ｃ＋＋代码）
#include
#include
intmain（）
{
doublea,b,n;
doubleh,S,x;
intk=1;
while（scanf（"%lf%lf%lf",&a,&b,&n））
{
h=（b-a）/n;
S=sin（b）/b-sin（a）/a;
x=a;
while（k<=n）
{
x=h/2+x;
S+=4*sin（x）/x;
x+=h/2;
S+=2*sin（x）/x;
k++;
}
S=S*h/6;
printf（"%lf\n",S）;
}
return0;
}
升级会员 