牛顿Word文档格式.docx

牛顿Word文档格式.docx

《牛顿Word文档格式.docx》由会员分享，可在线阅读，更多相关《牛顿Word文档格式.docx（11页珍藏版）》请在冰豆网上搜索。

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三个参数

%n，牛顿-科特斯数公式的阶数，必须满足1<

n<

7，因为n>

=8时不能保证公式的稳定性

（1）n=1，即梯形公式

（2）n=2，即辛普森公式

（3）n=4，即科特斯公式

%y，数值积分结果

%Ck，科特斯系数

%Ak，求积系数

fun1=@（x）sin（x）;

%必须可以接受矢量输入

%fun2=[0:

0.1:

0.5;

sin（0:

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

fun构成应为：

第一列为x，第二列为y，并且个数为小于10的等距节点！

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

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

输入参数错误，请参考函数帮助！

%计算科特斯系数

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阶科特斯系数

n+1

k=i-1;

Ck（i）=（-1）^（n-k）/factorial（k）/factorial（n-k）/n*quadl（@（t）intfun（t,n,k）,0,n）;

functionf=intfun（t,n,k）

%科特斯系数中的积分表达式

f=1;

fori=[0:

k-1,k+1:

n]

f=f.*（t-i）;

求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<

stdio.h>

math.h>

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++）{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;

n;

i++）{p2+=Cotes（a+i*h,a+（i+1）*h）;

returnp2;

intmain（）

doublearea;

area=Cotes（0,8）;

printf（"

Cotes:

area=%f\n"

area）;

area=Complex（0,8,0.0001）;

Complex:

area=Complex（0,8,0.00000001）;

return1;

C++编程用梯形求积公式求解定积分∫3lnxdx积分区间为（1，2）的值

#include"

stdio.h"

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;

}

%lf\n"

s）;

iostream>

cmath>

usingnamespacestd;

intmain（）{

doubleh=1e-6;

doublesum=0;

2）{

sum+=log（x）;

x+=h;

sum=sum*h*3;

cout<

<

sum<

endl;

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<

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;

S）;

return0;