实验二 定积分的近似计算

定积分近似计算 问题背景和实验目的 定积分计算的基本公式是牛顿－莱布尼兹公式。但当被积函数的原函数不知道时，如何计算？这时就需要利用近似计算。特别是在许多实际应用中，被积函数甚至没有解析表达式，而是一条实验记录曲线，或一组离散的采样值，此时只能用近似方法计算定积分。 本实验主要介绍计算定积分的三种基本近似算法：矩形法、梯形法和抛物线法。同时介绍 Matlab 计算定积分的相关函数。

主要内容 计算定积分的基本数值算法 Matlab 计算积分的相关函数 矩形法 梯形法 抛物线法 数值积分函数 trapz、quad、integral、integral2 符号积分函数：int

定积分的近似 定积分的定义 n 充分大，x 充分小

矩形法 取 步长 节点 点 的常见取法：左端点，右端点或中点 左 点 法 右 点 法 中 点 法

矩形法举例 例： 解： h =1/n=0.01, xi = i*h, a=0, b=1, n=100 fuluA.m 解： h =1/n=0.01, xi = i*h, a=0, b=1, n=100 (i = 0, 1, 2, ..., 100) 左点法： 右点法： 中点法：

补：函数句柄 函数句柄的定义 函数句柄：可以理解成一个函数的代号或别名，调用函数句柄就等价于调用该函数。 fhandle = @ 函数名 @ 的作用就是将一个函数的函数句柄赋值给左边的变量 例： f = @sin; y = f(pi/3)

补：匿名函数 匿名函数是 MATLAB的一种函数描述形式，可以让用户编写简单的函数而不需要创建M文件，并且具有很高的执行效率。 匿名函数的定义 f = @ (变量列表) 表达式 这里返回的 f 是一个函数句柄（function handle） 例： f = @(x) 1/(1+x*x); % 一个自变量 y = f(2) whos f % f 是一个函数句柄 f = @(x,y) x^2 + y^2; % 两个自变量 y = f(2,3) % 注意对应关系

补：匿名函数 匿名函数中可以采用数组运算，以便作用在向量或数组上 例： f = @(x) 1./(1+x.*x); % 一个自变量 y1 = f(x1); % f 是作用在向量 x1 的每个分量上 如果需要同时计算函数在多个点上的值，可以将函数直接作用在向量上，但此时函数的定义中必须采用数组运算！

矩形法举例 相对误差分析 不同的算法有不同的计算精度 理论值： 左点法相对误差： 右点法相对误差： 中点法相对误差： 有没有更好的近似计算定积分的方法 ?

定积分几何意义

梯形法 曲边小梯形的面积可以由直边小梯形的面积来近似 整个曲边梯形的面积：

梯形法 如果我们 n 等分区间 [a,b]，即令： 则 梯形公式 梯形公式与中点公式有什么区别 ?

梯形法举例 例：用梯形法计算定积分 ( 取 n=100 )，并计算相对误差 解： ==> fuluB.m 解： a=0, b=1, n=100, f (x) = 1/( 1+x2 ) ==> h =1/100=0.01, xi = i*h, yi = f (xi) ==> 相对误差：

抛物线法 n 等分区间 [a,b] ，得 计算节点和中点上的函数值： 在区间 [xi-1, xi] 上，用过以下三点 用抛物线代替该直线，计算精度是否会更好？ 的抛物线来近似原函数 f (x) 。

抛物线法 设过以上三点的抛物线方程为： y =  x2 +  x +  = pi (x) 则在区间 [xi-1, xi] 上，有

抛物线法 相加后可得： 抛物线法公式 或 辛卜生 (Simpson) 公式

抛物线法 例：用抛物线法计算下面定积分 ( 取 n=100 )，并计算相对误差 解： a=0, b=1, n=100, 相对误差： fuluC.m 解： a=0, b=1, n=100, 相对误差：

主要内容 计算定积分的基本数值算法 Matlab 计算积分的相关函数 矩形法 梯形法 抛物线法 数值积分函数 trapz、quad、integral、integral2 符号积分函数：int

trapz 梯形法 trapz(x, y) x 为分割点（节点）组成的向量， y 为被积函数在节点上的函数值组成的向量。

trapz 举例 例：用梯形法计算下面定积分 ( 取 n=100 ) 解： a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 ) 前面的做法 trapz函数 x=0:1/100:1; y=1./(1+x.^2); inum=trapz(x, y)

quad 自适应抛物线法 quad(f,a,b,tol) f 是函数句柄，也可用字符串表示（不推荐）， 其中涉及的运算必须采用数组运算 f = f(x) 为被积函数，[a,b] 为积分区间，tol 为计算精度 不用自己分割积分区间 可以指定计算精度，若不指定，缺省精度是 10-6 精度越高，函数运行的时间越长 f 是函数句柄，也可用字符串表示（不推荐）， 其中涉及的运算必须采用数组运算 将自变量看成是向量！

quad 举例 例：用 quad 计算定积分： 解： f=@(x) 1./(1+x.^2); inum=quad(f, 0, 1) % 采用缺省精度 inum=quad(@(x) 1./(1+x.^2), 0, 1, 1e-10)

integral 全局自适应积分法（R2012a以后版本） integral(f,a,b) integral(f,a,b,'RelTol',tol) 该函数比 quad 效率更高，且可以处理一些非正常积分 可以指定计算精度，若不指定，缺省精度是 10-6 f 必须是函数句柄，且涉及的运算必须采用数组运算 f=@(x) 1./(1+x.^2); inum=integral(f,0,1) inum=integral(f,0,1,'RelTol',1e-10) f=@(x) exp(-x); inum=integral(f,0,inf)

integral2 计算二重积分的全局自适应积分法 integral2(f,a,b,c,d,tol) integral2(f,a,b,c,d,'RelTol',tol) 可以指定计算精度，若不指定，缺省精度是 10-6 f 必须是函数句柄，且涉及的运算必须采用数组运算

integral2 例：计算二重积分 f=@(x,y) 4*x.*y+3*y.^2; inum=integral2(f,-1,1,0,2) 在前面的是第一积分变量，在后面的是第二积分变量 f=@(x,y) 4*x.*y+3*y.^2; inum=integral2(f,-1,1,0,2) 注意积分变量与积分区间的对应关系

int 符号积分 int(f,v,a,b) % 计算定积分 int(f,a,b) % 计算关于默认变量的定积分 syms x; f=1/(1+x^2); inum=int(f,x,0,1)

相关函数 double(a) 将 a 转化为双精度型，若 a 是字符，则取对应的 ASCII 码 例： a=3; x=3 x=double(a) y=double('a') x=3 y=97

数值实验 例：用 Matlab 函数近似计算定积分 x=1:0.001:2; y=exp(x.^(-2)); inum=trapz(x,y) 梯形法： x=1:0.001:2; y=exp(x.^(-2)); inum=trapz(x,y) 抛物线法： f=@(x) exp(x.^(-2)); inum=quad(f, 1, 2, 1e-10) 用Matlab演示 符号积分法： syms x; inum=int(exp(x^(-2)),x,1,2)

数值实验 例：用 Matlab 函数近似计算二重积分 f=@(x,y) x+y.^2; 数值积分法： f=@(x,y) x+y.^2; inum=integral2(f, 0, 2, -1, 1) 符号积分法： 用Matlab演示 syms x y; f=int(x+y^2,y,-1,1); inum=int(f,x,0,2)

上机作业 上机作业 要求填写实验报告 要求： 没有提到使用Matlab函数的题必须使用编程完成 教材第 75 页：2、3、6 将所编写的程序分别命名为 hw221.m, hw222.m, hw231.m, hw261.m, hw262.m 要求： 用Matlab演示 没有提到使用Matlab函数的题必须使用编程完成 要求使用Matlab函数或命令的必须使用函数或命令 使用课程主页上的附录程序

上机要求 上机要求 将所有文件作为附件，通过 foxmail 以邮件形式发给 mhjs@system.mail 邮件主题为：机号-学号-姓名，其中机号为 两位数 三个字段之间用英文状态下的减号连接 每个 M 文件的第一行添加一条注解语句: % 机号-学号-姓名