实验二定积分的近似计算

利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．

本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，Matlab有专门函数可用．

1．sum(a)：求数组a的和．

2．format long：长格式，即屏幕显示15位有效数字．

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

3．double()：若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值．

4．quad()：抛物线法求数值积分．

格式：quad(fun,a,b)，注意此处的fun是函数，并且为数值形式的，所以使用*、/、^等运算时要在其前加上小数点，即.*、./、.^等．

例：Q = quad('1./(x.^3-2*x-5)',0,2);

5．trapz()：梯形法求数值积分．

格式：trapz(x,y)

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

例：计算

x=0:pi/100:pi;y=sin(x);

trapz(x,y)

6．dblquad()：抛物线法求二重数值积分．

格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun可以用inline定义，也可以通过某个函数文件的句柄传递．

例1：Q1 = dblquad(inline('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：

function z = integrnd(x, y)

z = y*sin(x);

7．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．syms变量1变量2…：定义变量为符号．

9．sym('表达式')：将表达式定义为符号．

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

10．int(f,v,a,b)：求f关于v积分，积分区间由a到b．

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

1．矩形法

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

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

针对不同的取法，计算结果会有不同，我们以为例(取)，

(1)左点法：对等分区间

，

在区间上取左端点，即取，

0.78789399673078，

理论值，此时计算的相对误差

(2)右点法：同(1)中划分区间，在区间上取右端点，即取，

0.78289399673078，

理论值，此时计算的相对误差

(3)中点法：同(1)中划分区间，在区间上取中点，即取，

0.78540024673078，

理论值，此时计算的相对误差

如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物线法就是这一指导思想的产物．

2．梯形法

等分区间

，

相应函数值为

()．

曲线上相应的点为

()

将曲线的每一段弧用过点，的弦(线性函数)来代替，这使得每个上的曲边梯形成为真正的梯形，其面积为

，．

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

，

即

，

称此式为梯形公式．

仍用的近似计算为例，取，

0.78539399673078，

理论值，此时计算的相对误差

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

3．抛物线法

由梯形法求近似值，当为凹曲线时，它就偏小；当为凸曲线时，它就偏大．若每段改用与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法．

将积分区间作等分，分点依次为

，，

对应函数值为

()，

曲线上相应点为

()．

现把区间上的曲线段用通过三点，，的抛物线

来近似代替，然后求函数从到的定积分：

由于，代入上式整理后得

同样也有

……

将这个积分相加即得原来所要计算的定积分的近似值：

，

即

这就是抛物线法公式，也称为辛卜生(Simpson)公式．

仍用的近似计算为例，取，

=0.78539816339745，

理论值，此时计算的相对误差

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

(1)数值计算

方法1：int('1/(1+x^2)','x',0,1)    (符号求积分)

方法2：quad('1./(1+x.^2)',0,1)(抛物线法求数值积分)

方法3：x=0:0.001:1;

y=1./(1+x.^2);

trapz(x,y)(梯形法求数值积分)

(2)数值计算

方法1：int(int('x+y^2','y',-1,1),'x',0,2)(符号求积分)

方法2：dblquad(inline('x+y^2'),0,2,-1,1)(抛物线法二重数值积分)

1．实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算，取，并比较三种方法的精确程度．

2．分别用梯形法与抛物线法，计算，取．并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异．

3．试计算定积分．(注意：可以运用trapz()、quad()或附录程序求解吗？为什么？)

4．将的近似计算结果与Matlab中各命令的计算结果相比较，试猜测Matlab中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点．

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

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