实验二

实验二 定积分的近似计算

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

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

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

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

解： a=0, b=1, n=100 h =1/n=0.01, xi = i*h, (i = 0, 1, 2, ..., 100) • 左点法： • 右点法： • 中点法：矩形法举例例：用三种不同的矩形法计算下面的定积分 ( 取 n=100 )，并比较这三种方法的相对误差。 fuluA.m

补：函数句柄 函数句柄：可以理解成一个函数的代号或别名，调用函数句柄就等价于调用该函数。 • 函数句柄的定义 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); % 一个自变量 x1 = 0 : 0.01 : 1; % x1 是向量 y1 = f(x1); % f 是作用在向量 x1 的每个分量上如果需要同时计算函数在多个点上的值，可以将函数直接作用在向量上，但此时函数的定义中必须采用数组运算！

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

定积分几何意义

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

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

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

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

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

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

相对误差： 抛物线法例：用抛物线法计算下面定积分 ( 取 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(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) 注意积分变量与积分区间的对应关系

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

相关函数 double(a) 将a转化为双精度型，若a是字符，则取对应的ASCII 码例： a=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) • 抛物线法： f=@(x) exp(x.^(-2)); inum=quad(f, 1, 2, 1e-10) • 符号积分法： syms x; inum=int(exp(x^(-2)),x,1,2)

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

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

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

实验二

实验二

Presentation Transcript