定积分的近似计算

1、实验二  定积分的近似计算一、问题背景与实验目的利用牛顿莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法对于定积分的近似数值计算，Matlab有专门函数可用 二、相关函数（命令）及简介1sum(a)：求数组a的和2format long：长格式，即屏幕显示15位有效

2、数字（注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）3double()：若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值4quad()：抛物线法求数值积分格式： quad(fun,a,b) ，注意此处的fun是函数，并且为数值形式的，所以使用*、/、等运算时要在其前加上小数点，即 .*、./、.等例：Q = quad('1./(x.3-2*x-5)',0,2);5trapz()：梯形法求数值积分格式：trapz(x,y)其中x为带有步长的积分区间；y为数值形式的运算（相当于上面介绍的函数fun）例：计算x=0:pi/100:p

3、i;y=sin(x);trapz(x,y)6dblquad()：抛物线法求二重数值积分格式：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)'

4、), 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);7fprintf（文件地址，格式，写入的变量）：把数据写入指定文件例：x = 0:.1:1;y = x; exp(x);fid = fopen('exp.txt','w');   %打开文件fprintf(fid,'%6.2f %12.8fn',y);  

5、60; %写入fclose(fid)                 %关闭文件8syms 变量1 变量2 ：定义变量为符号9sym('表达式')：将表达式定义为符号解释：Matlab中的符号运算事实上是借用了Maple的软件包，所以当在Matlab中要对符号进行运算时，必须先把要用到的变量定义为符号10int(f,v,a,b)：求f关于v积分，积分区间由a到b11subs(f，'x'，a)：将 a 的值

6、赋给符号表达式 f 中的 x，并计算出值若简单地使用subs(f)，则将f的所有符号变量用可能的数值代入，并计算出值。三、实验内容一、问题的提出 计算定积分的方法： (1) 求原函数； (2) 利用牛顿莱布尼茨公式计算结果。 问题： (1) 被积函数的原函数不能用初等函数表示； (2) 被积函数难于用公式表示，而是用图形或表格给出的； (3) 被积函数虽然能用公式表示，但计算其原函数很困难。 解决办法： 建立定积分的近似计算方法。 思路： 在数值上表示曲边梯形的面积，只要近似地算出相应的曲边梯形的面积，就可得到所给定积分的近似值。 常用方法：矩形法、梯形法、抛物线法 二、矩形法 用分

7、点a=x0,x1,xn=b将区间a,bn等分，取小区间左端点的函数值yi(i=0,1,2,n-1)作为窄矩形的高，如图： 则有： 取小区间右端点的函数值yi(i=1,2,n)作为窄矩形的高，如图： 则有： 以上两公式称为矩形法公式。 例：用矩形法求 ，并与用牛顿莱布尼茨公式计算的结果进行比较。 程序如下： #include <stdio.h> #include <math.h>  void main()     double result,a=0,b=1,i,n=1000000,h;

8、60;        printf("按牛顿公式计算得到的结果：%fn",sin(b)-sin(a);     result=0;     h=(b-a)/n;/计算区间高度     for(i=1;i<=n;i+)/求和         result=result+cos(a+i*h);  

9、   result=h*result;/乘以区间高度     printf("用近似公式计算得到的结果：%fn",result); 三、梯形法 梯形法就是在每个小区间上，以窄梯形的面积近似代替窄曲边梯形的面积，如图： 则有：例：用梯形法求 ，并与用牛顿莱布尼茨公式计算的结果进行比较。 #include <stdio.h> #include <math.h>  void main()     double result,a=

10、0,b=1,i,n=1000000,h;          printf("按牛顿公式计算得到的结果：%fn",sin(b)-sin(a);     result=0;     h=(b-a)/n;/计算区间高度     for(i=1;i<=n-1;i+)/求和         result

11、=result+cos(a+i*h);     result+=(cos(a)+cos(b)/2;     result=h*result;/乘以区间高度     printf("用近似公式计算得到的结果：%fn",result); 四、抛物线法 此法就是将曲线分成许多小段，用对称轴平行于y轴的二次抛物级上的一段弧来近似替代原来的曲线弧，从而得到定积分的近似值。 用分点a=x0,x1,xn=b将区间a,bn等分（偶数），这些分点对应曲线上的点为Mi(xi,

12、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

13、,h上，计算所得的定积分的值与原区间上的相同。 计算在-h,h上过三点的抛物线为曲边的面积。 抛物线中的可由下列方程组确定： 由此得：于是所求面积为： 显然，曲边梯形的面积只与的纵坐标及底边所在的区间的长度2h有关。由此可知n/2组梯形的面积为： 例：用抛物法求 ，并与用牛顿莱布尼茨公式计算的结果进行比较。 #include <stdio.h> #include <math.h>  void main()     double result,a=0,b=1,i,n=1000000,h;   

14、       printf("按牛顿公式计算得到的结果：%fn",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

15、*i-1)*h);     result+=cos(a)+cos(b);     result=h*result/3;/乘以区间高度     printf("用近似公式计算得到的结果：%fn",result); 注意：对于以上三种方法当n取得越大时近似程度就越好。 练习题： 用三种积分法近似计算如下定积分的值： 4. 直接应用Matlab命令计算结果 （1）    数值计算方法1：int('1/(1

16、+x2)','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+y2','y',-1,1),'x',0,2)   （符号求积分）方法2：dblquad(inline('x+y2'),0,2,-1,1)   （抛物线法二重数值积分） 四、自己动手1  实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算，取，并比较三种方法的精确程度2  分别用梯形法与抛物线法，计算，取并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异3