(整理)定积分的近似计算.

实验二定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，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) 求原函数；(2) 利用牛顿－莱布尼茨公式计算结果。

问题：(1) 被积函数的原函数不能用初等函数表示；(2) 被积函数难于用公式表示，而是用图形或表格给出的；(3) 被积函数虽然能用公式表示，但计算其原函数很困难。

解决办法：建立定积分的近似计算方法。

思路：在数值上表示曲边梯形的面积，只要近似地算出相应的曲边梯形的面积，就可得到所给定积分的近似值。

常用方法：矩形法、梯形法、抛物线法．二、矩形法用分点a=x0,x1, (x)n=b将区间[a,b]n等分，取小区间左端点的函数值yi(i=0,1,2,…,n-1)作为窄矩形的高，如图：则有：取小区间右端点的函数值y(i=1,2,…,n)作为窄矩形的高，如图：i则有：以上两公式称为矩形法公式。

例：用矩形法求，并与用牛顿－莱布尼茨公式计算的结果进行比较。

程序如下：#include <stdio.h>#include <math.h>void main(){double result,a=0,b=1,i,n=1000000,h;printf("按牛顿公式计算得到的结果：%f\n",sin(b)-sin(a));result=0;h=(b-a)/n;//计算区间高度for(i=1;i<=n;i++)//求和result=result+cos(a+i*h);result=h*result;//乘以区间高度printf("用近似公式计算得到的结果：%f\n",result);}三、梯形法梯形法就是在每个小区间上，以窄梯形的面积近似代替窄曲边梯形的面积，如图：则有：例：用梯形法求，并与用牛顿－莱布尼茨公式计算的结果进行比较。

#include <stdio.h>#include <math.h>void main(){double result,a=0,b=1,i,n=1000000,h;printf("按牛顿公式计算得到的结果：%f\n",sin(b)-sin(a));result=0;h=(b-a)/n;//计算区间高度for(i=1;i<=n-1;i++)//求和result=result+cos(a+i*h);result+=(cos(a)+cos(b))/2;result=h*result;//乘以区间高度printf("用近似公式计算得到的结果：%f\n",result);}四、抛物线法此法就是将曲线分成许多小段，用对称轴平行于y轴的二次抛物级上的一段弧来近似替代原来的曲线弧，从而得到定积分的近似值。

用分点a=x0,x1, (x)n=b将区间[a,b]n等分（偶数），这些分点对应曲线上的点为M i (xi,yi)(其中yi=f(xi)，i=0,1,2,…,n)，如图：因为经过三个不同的点可以唯一确定一条抛物线，可将这些曲线上的点Mi互相衔接地分成n/2组{M0,M1,M2},{M2,M3,M4},…,{Mn-2,Mn-1,Mn}，即每相邻两个区间为一组。

在每组{ M2k-2, M2k-1, M2k}(k=1,2,…,n/2)所对应的子曲间[x2k-2,x2k]上，用经过点M2k-2, M2k-1, M2k的二次抛物线近似代替曲线弧。

下面讨论如何计算积分。

设h为区间高度，即h=x2k -x2k-1=x2k-1-x2k-2。

根据积分性质（积分在数值上表示曲边梯形的面积）有如下等式成立：即将区间[x2k-2,x2k]平移到区间[-h,h]上，计算所得的定积分的值与原区间上的相同。

计算在[-h,h]上过三点的抛物线为曲边的面积。

抛物线中的可由下列方程组确定：由此得：于是所求面积为：显然，曲边梯形的面积只与的纵坐标及底边所在的区间的长度2h 有关。

由此可知n/2组梯形的面积为：例：用抛物法求，并与用牛顿－莱布尼茨公式计算的结果进行比较。

#include <stdio.h>#include <math.h>void main(){double result,a=0,b=1,i,n=1000000,h;printf("按牛顿公式计算得到的结果：%f\n",sin(b)-sin(a));result=0;h=(b-a)/n;//计算区间高度for(i=1;i<=n/2;i++)//求和result=result+2*cos(a+2*i*h)+4*cos(a+(2*i-1)*h);result+=cos(a)+cos(b);result=h*result/3;//乘以区间高度printf("用近似公式计算得到的结果：%f\n",result);}注意：对于以上三种方法当n取得越大时近似程度就越好。

练习题：用三种积分法近似计算如下定积分的值：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 循环。