牛顿.docx

1、牛顿牛顿-柯特斯数值积分公式及其MATLAB的实现(2009-11-20 23:46:18)至于Newton-Cotes数值积分的具体表达、推导、说明和附件，下面给只是给出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

2、-Cotes公式的阶数，必须满足n8，否则积分没法保证稳定性%(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);for i=1:ms(i)=NewtonCotes(fun,xk(i),xk(i+1),n);endy=sum(s);牛顿-科特斯数值积分公式function y,Ck,Ak=NewtonCotes(fun,a,b,n)% y=NewtonCotes(fun,a,b,n)% 牛顿-

3、科特斯数值积分公式% 参数说明：% 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，牛顿-科特斯数公式的阶数，必须满足1n=8时不能保证公式的稳定性%(1)n=1，即梯形公式%(2)n=2，即辛普森公式%(3)n=4，即科特斯公式% y，数值积分结果% Ck，科特斯系数

4、% 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% contact mematlabsky% 2009-11-20 15:06:51%if nargin=1mm,nn=size(fun);if mm=8error(为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！)elseif

5、 nn=2error(fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！)endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm-1;elseif nargin=4% 计算积分节点xk和节点函数值fxxk=linspace(a,b,n+1);if isa(fun,function_handle)fx=fun(xk);elseerror(fun积分函数的句柄，且必须能够接受矢量输入！)endelseerror(输入参数错误，请参考函数帮助！)end% 计算科特斯系数Ck=cotescoeff(n);% 计算求积系数Ak=(b-a

6、)*Ck;% 求和算积分y=Ak*fx;function Ck=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阶科特斯系数for i=1:n+1k=i-1;Ck(i)=(-1)(n-k)/

7、factorial(k)/factorial(n-k)/n*quadl(t)intfun(t,n,k),0,n);endfunction f=intfun(t,n,k)% 科特斯系数中的积分表达式f=1;for i=0:k-1,k+1:nf=f.*(t-i);end求f(x)在a,b上的定积分用柯特斯公式把a,b 4等分 步长h=(b-a)/4 取等分点xi=a+i*h (i=0,1,2,3,4)柯特斯公式为C=1/90 *(b-a)*(7*f(x0)+32*f(x1)+12*f(x2)+32*f(x3)+7*f(x4)复化求积分先将a,b 等分成n份 步长为h=(b-a)/n ,xk=a+k

8、*h (k=0,1,2.n)先用柯特斯公式在每个子段xk,xk+1 求得积分 ck然后求数列 c0,c1,-1的和作为积分 zhulinpptor (zhulin) 等级： #include #include typedef double (*pFun)(double ); double fun(double x)return x*x*x*x*x*x;double Cotes(double a,double b)/柯特斯公式 有5次代数精度/如果机械求积公式对xk均能准确成立，/但对不准确x(k+1)则称机械求积公式具有k次代数精度double c5;double value=0.0;int

9、d5=7,32,12,32,7;double h=(b-a)/4.0;for(int i=0;i5;i+)ci=a+i*h;pFun f=fun;for(i=0;ie) p1=p2; n=2*n;h=(b-a)/n;p2=0;for(i=0;in;i+)p2+=Cotes(a+i*h,a+(i+1)*h); return p2;int main() double area; area=Cotes(0,8); printf(Cotes:area=%fn,area); area=Complex(0,8,0.0001); printf(Complex:area=%fn,area); area=Com

10、plex(0,8,0.00000001); printf(Complex:area=%fn,area); return 1;C+编程 用梯形求积公式求解定积分3lnxdx积分区间为（1，2）的值 #include stdio.h#include math.h#define N 100void main() double delta=1.0/N; double x=1; double y1,y2; double s=0; y1=0; while(x2) y2=3*log(x+delta); s+=(y1+y2)*delta/2; y1=y2; x+=delta; printf(%lfn,s);#

11、include #include using namespace std;int main() double x=1; double h=1e-6; double sum=0; while(x2) sum+=log(x); x+=h; sum=sum*h*3; coutsum chazhijifenx=127.5040tixing 积分精确值Iexact=3.141592653 n I Error 2 3.100000000 0.041592653 4 3.131176471 0.010416182 8 3.138988494 0.002604159 16 3.140941612 0.0006

12、51041 32 3.141429893 0.000162760 64 3.141551963 0.000040690128 3.141582481 0.000010172256 3.141590110 0.000002543512 3.141592018 0.000000635simpson Simpson法计算积分 精确值Iexact=3.141592653;n I Error 2 3.133333333 0.008259320 4 3.141568627 0.000024026 8 3.141592502 0.000000151 16 3.141592651 0.000000002 32

13、 3.141592654 -0.000000001 64 3.141592654 -0.000000001128 3.141592654 -0.000000001256 3.141592654 -0.000000001512 3.141592654 -0.000000001复化辛甫生求积公式（代码）#include#includeint main()double a,b,n;double h,S,x;int k=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(%lfn,S);return 0;