定积分的近似计算.docx

定积分的近似计算.docx

《定积分的近似计算.docx》由会员分享，可在线阅读，更多相关《定积分的近似计算.docx（11页珍藏版）》请在冰点文库上搜索。

定积分的近似计算

实验

定积分的近似计算

一、问题背景与实验目的

利用牛顿一莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被

积函数的原函数能用初等函数表达出来的情形.如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法.在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这

时只能应用近似方法去计算相应的定积分.

本实验将主要研究定积分的三种近似计算算法：

矩形法、梯形法、抛物线法.对于定积分的近似数值计算，Matlab有专门函数可用.

二、相关函数（命令）及简介

1.sum（a）：

求数组a的和.

2.formatlong：

长格式，即屏幕显示15位有效数字.

（注：

由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）

3.doubleO：

若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值.

4.quad（）：

抛物线法求数值积分.

格式：

quad（fun,a,b），注意此处的fun是函数，并且为数值形式的，所以使用*、人八等运算时要在其前加上小数点，即.*、./、A等.

例：

Q=quad（'1./（xA3-2*x-5）',0,2）;

5.trapz（）：

梯形法求数值积分.

格式：

trapz（x,y）

其中x为带有步长的积分区间；y为数值形式的运算（相当于上面介绍的函数fun）

例：

计算『sin（x）dx

x=0:

pi/100:

pi;y=sin（x）;

trapz（x,y）

6.dblquad（）：

抛物线法求二重数值积分.

格式：

dblquad（fun,xmin,xmax,ymin,ymax），fun可以用inline定义，也可以通过某个函数文件的句柄传递.

例1：

Q1=dblquad（iniine（'y*sin（x）'）,pi,2*pi,0,pi）

顺便计算下面的Q2,通过计算，比较Q1与Q2结果（或加上手工验算），找出积分变量x、y的上下限的函数代入方法.

Q2=dblquad（inline（'y*sin（x）'）,0,pi,pi,2*pi）

例2：

Q3=dblquad（@integrnd,pi,2*pi,0,pi）

这时必须存在一个函数文件integrnd.m：

functionz=integrnd（x,y）

7.

z=y*sin（x）;

fprintf（文件地址，格式，写入的变量）：

把数据写入指定文件.

例:

x=0:

.1:

1;

y=[x;exp（X）];

fid=fopen（'exp.txt','w'）;%打开文件

fprintf（fid,'%6.2f%12.8f\n',y）;%写入

fclose（fid）%关闭文件

8.

9.

syms变量1变量2…：

定义变量为符号.

sym（'表达式'）：

将表达式定义为符且

解释：

Matlab中的符号运算事实上是借用了Mapie的软件包，所以当在Matlab中要对符号进行运算时，必须先把要用到的变量定义为符号.

10.int（f,v,a,b）:

求f关于v积分，积分区间由a到b.

11.subs（f,'x',a）：

将a的值赋给符号表达式f中的x，并计算出值.若简单地使用subs（f），则将f的所有符号变量用可能的数值代入，并计算出值.

三、实验内容

1.矩形法

根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即

bn严

Lf（x）dx=Sf（UjAxi

ay

在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法.不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度.

针对不同\的取法，计算结果会有不同，我们以r伞为例（取n=100），

01+x

（1）左点法：

对等分区间

X0=aCX1V…

n

在区间[Xi4，Xi]上取左端点，即取E-Xy,

1dxn

——=Zf

（2）也Xi止0.78789399673078

1+x匕

（2）右点法：

同

（1）中划分区间，在区间［Xi_L,Xi］上取右端点，即取H=Xi，

1dxn

——=2f（勺疋x—0.78289399673078

1+Xy

理论值J0£r4'此时计算的相对误差

0.7854002467307〜兀/4

止2.653x10厘

如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那

么可以期望得到比矩形法效果好得多的近似计算公式.下面介绍的梯形法和抛物

线法就是这一指导思想的产物.

2.梯形法

等分区间

b—a.b—a

Xo=a<%<…

nn

相应函数值为

y0,yi,…,yn（y^f（Xi）,i=0,1，…,n）.

曲线y=f（X）上相应的点为

P0,R,…，Pn（P=（Xi,yi）,i=0,1，…，n）

将曲线的每一段弧P^Pj用过点Rd，P的弦P/R（线性函数）来代替，这使得每个［Xi」，Xi］上的曲边梯形成为真正的梯形，其面积为

y』+yiX虫X，i=12…，n.

2

于是各个小梯形面积之和就是曲边梯形面积的近似值,

bny+y氐Xn

Jaf（x）dx止2y*|咒3=〒无（yi_L+yi），

i422y

Jabf（x）dx止牛（号+yi+卅+y」+普），

称此式为梯形公式.

1dx

仍用f竺y的近似计算为例，取n=100，

，01+x2

dX^-_a（匹+%+山+yn」中一）=0.78539399673078,1+xn22

理论值［1-dx^=-，此时计算的相对误差

‘01+x24

很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多.

3.抛物线法

由梯形法求近似值，当y=f（x）为凹曲线时，它就偏小；当y=f（x）为凸曲

线时，它就偏大.若每段改用与它凸性相接近的抛物线来近似时，就可减少上述

缺点，这就是抛物线法.

将积分区间［a,b］作2n等分，分点依次为

x0=acxrV

2n2n

对应函数值为

yo,yi,…，y2n（y^f（Xi）,i=o,i，…,2n）,

曲线上相应点为

Po,R,…，P2n（P=（Xi,yi）,i=0,1，…，2n）.

P2（X2,y2）的抛物线

y=ax2+Px+Y=5（x）

来近似代替，然后求函数Pi（X）从Xo到X2的定积分:

x?

X22Ot33p22

fPi（x）dx=J（otx2+Px+Y）dx=—（X2-Xo）+—（X2-x0）+丫化—x0）

・Xo“32

XX

[（ax0+Px0+Y）+（ax；+Px?

+丫）+a（x0+x2）2+2P（x0+X2）+4丫]

X2

尹（yo+4y1+y2）=詈（yo+4y1+y2）

同样也有

x4b—a

Jx2P2（X）dX=^（y2+4y3+y4）

1dxb-a

贡[y。

+y2n+4（y1+y3川+沁）+2（y2+y4仙+皿）]

=0.78539816339745,

理论值［1—吗=兰，此时计算的相对误差

‘01+x24

4.直接应用Matlab命令计算结果

四、自己动手

1dx

1.

实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算f羊,

、01+X

取n=258，并比较三种方法的精确程度.

2dx

2.

分别用梯形法与抛物线法，计算fdx，取n=120.并尝试直接使用函数

‘1x

trapz（）、quad（）进行计算求解，比较结果的差异.

3.

试计算定积分f讫沁dx.（注意：

可以运用trapz（）、quad（）或附录程序求解

■0x

吗？

为什么？

）

1dx

4.

将J二的近似计算结果与Matlab中各命令的计算结果相比较，试猜测

01+x

Matlab中的数值积分命令最可能采用了哪一种近似计算方法？

并找出其他例子支持你的观点.

通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型的函数（具有某种单调特性或凹凸特性），用某种近似计算方法所得结果更接

近于实际值？

6.学习fulu2sum.m的程序设计方法，尝试用函数sum改写附录1和附录3的

程序，避免for循环．

五、附录

附录1：

矩形法（左点法、右点法、中点法）（fulu1.m）formatlongn=100;a=0;b=1;

inum1=0;inum2=0;inum3=0;symsxfx

fx=1/（1+x^2）;

%左点

%右点

%左点值%右点值%中点值

fori=1:

nxj=a+（i-1）*（b-a）/n;xi=a+i*（b-a）/n;

fxj=subs（fx,'x',xj）;fxi=subs（fx,'x',xi）;fxij=subs（fx,'x',（xi+xj）/2）;inum1=inum1+fxj*（b-a）/n;inum2=inum2+fxi*（b-a）/n;inum3=inum3+fxij*（b-a）/n;

endinum1inum2inum3integrate=int（fx,0,1）integrate=double（integrate）fprintf（'Therelativeerrorbetweeninum1andreal-valueisabout:

%d\n\n',...

abs（（inum1-integrate）/integrate））fprintf（'Therelativeerrorbetweeninum2andreal-valueisabout:

%d\n\n',...

abs（（inum2-integrate）/integrate））

附录2：

梯形法（fulu2.m）formatlongn=100;a=0;b=1;inum=0;

symsxfxfx=1/（1+x^2）;

fori=1:

nxj=a+（i-1）*（b-a）/n;xi=a+i*（b-a）/n;fxj=subs（fx,'x',xj）;fxi=subs（fx,'x',xi）;inum=inum+（fxj+fxi）*（b-a）/（2*n）;

end

inum

integrate=int（fx,0,1）integrate=double（integrate）fprintf（'Therelativeerrorbetweeninumandreal-valueisabout:

%d\n\n',...

abs（（inum-integrate）/integrate））

附录2sum：

梯形法（fulu2sum.m）,利用求和函数，避免for循环

formatlong

n=100;a=0;b=1;

symsxfx

fx=1/（1+x^2）;

%所有左点的数组

%所有右点的数组

%所有左点值

%所有右点值

%梯形面积%加和梯形面积求解

i=1:

n;

xj=a+（i-1）*（b-a）/n;

xi=a+i*（b-a）/n;

fxj=subs（fx,'x',xj）;

fxi=subs（fx,'x',xi）;

f=（fxi+fxj）/2*（b-a）/n;inum=sum（f）integrate=int（fx,0,1）integrate=double（integrate）

附录3：

抛物线法（fulu3.m）formatlongn=100;a=0;b=1;inum=0;

symsxfxfx=1/（1+x^2）;

%左点

%右点

%中点

fori=1:

n

xj=a+（i-1）*（b-a）/n;xi=a+i*（b-a）/n;

xk=（xi+xj）/2;fxj=subs（fx,'x',xj）;fxi=subs（fx,'x',xi）;fxk=subs（fx,'x',xk）;

inum=inum+（fxj+4*fxk+fxi）*（b-a）/（6*n）;endinumintegrate=int（fx,0,1）integrate=double（integrate）fprintf（'Therelativeerrorbetweeninumandreal-valueisabout:

%d\n\n',...