2 2 2 Ξ 2008 年 12 月 第 20 卷第 6 期 文章编号 :1009 4873 (2008) 06 石家庄职业技术学院学报 Journal of Shijiazhuang Vocational Technology Institute Dec. 2008 Vol. 20 　No. 6 0058 03 用 MATLAB 求数值积分的方法 陈佩宁a , 　刘 　竞b (石家庄职业技术学院 a. 信息工程系 ; b. 机电工程系 ,河北 石家庄 　050081) 摘 　要 :介绍了数值积分法的几种计算公式及相应的 MA TLAB 命令 ,并给出了用 MA TLAB 编程求数值积分 的实例. 关键词 :MA TLAB ;数值积分 ;矩形公式 ;梯形公式 ;辛普森公式 中图分类号 :O172 　　　文献标识码 :A 1 　引言 在一元微积分学中 , 若已知函数 f ( x) 在闭区 间[ a , b ] 上连续且其原函数为 F ( x) ,求 f ( x) 在该 区间上的定积分可用牛顿 - 莱布尼兹公式求解 , 即 ∫b MA TLAB 中可以用符号积分命令 int 求∫b 该命令格式为 : - F ( a) . 而 在 f ( x ) d x = F ( x ) | b a = F ( b) f ( x ) d x , a a int ( f , x , a , b) 　% 求函数 f 在区间[ a , b ] 上的 定积分. π 2 0 si nx dx. 例 1 　求∫ 解 　输入命令 : > > sy ms x > > I = i nt ( si n (x) ,x ,0 ,pi/ 2) , 结果显示为 : I = 1. 用牛顿 - 莱布尼兹公式计算定积分的方法在 理论上和解决实际问题中起到了很大的作用 ,但它 并不能解决定积分计算的所有问题. 在工程技术领 域常遇到十分复杂的情况而无法用牛顿 - 莱布尼 兹公式求解. 其可能出现的情况[ 1 ] 有 : (1) 某些被积函数 f (x) ,其原函数无法用初等 函数表示 ,如∫ex 2 dx ,∫si nx x dx 等. (2) 函数 f (x) 结构复杂 ,求其原函数非常困难. (3) 函数 f (x) 的结构虽然简单且其原函数存 (4) 函数 f (x) 没有具体的表达式 ,只有一些由 试验测试数据形成的表格或图形. 而在这些情况下 ,可采用“数值积分”的方法求 出定积分 (近似值) . 2 　数值积分法 用数值积分的方法求一个函数在区间[ a ,b ] 上 n n a , k = 1 n →∞ ∑ f (x) dx = li m 的定积分 ,可利用定积分的定义来求解 : f (ξk) b - a n I =∫b 设 In = ∑ 此时称 In 为数值积分. 显然 ,数值积分 In 就是 I 的近似值 ,并且当 n 越大 , In 就越接近于精确值 I. 由 于ξk 取值不同 ,数值积分 In 的结果会有所不同. [ 2 ] 数值积分的计算公式也有多种 : f (ξk) b - a n ,则 I = li m n →∞ In. k = 1 (1) 矩形公式 将积分区间[ a ,b ]n 等分 ,每个小区间宽度均为 h = (b - a) / n ,h 称为积分步长. 记 a = x0 < x1 < … < xk … < xn = b ,在小区 间上用小矩形面积近似小曲边梯形的面积 ,若分别 取左端点和右端点的函数值为小矩形的高 ,则分别 得到两个曲边梯形面积的近似计算公式 : L n = h ∑ k = 0 n- 1 f (xk) ,h = b - a n ① n f (xk) ,h = Rn = h ∑ 称公式 ①, ②分别为左 、右矩形公式 ,两个矩形 k = 1 n ② b - a 在 ,但其原函数的结构相对复杂. 面积分别小于和大于所求曲边梯形的面积. 收稿日期 :2007 作者简介 :陈佩宁 (1971 24 12 ) ,女 ,河北望都人 ,石家庄职业技术学院讲师. 

第 6 期 陈佩宁等 :用 MA TL A B 求数值积分的方法 95 (2) 梯形公式 如果将二者求平均值 ,则每个小区间上的小矩 形变为小梯形 ,整个区间上的值变为 : Tn = h ∑ f (xk) + n- 1 k = 1 h 2 [f (x0) + f (xn) ] ③ b - a 　　其中 ,h = . ③式称为梯形公式. n (3) 辛普森公式 为了提高计算结果的精度 ,用分段二次插值函 数代替 f (x) . 由于每段要用到相邻两个小区间端点 的三个函数值 ,所以 ,小区间的数目必须是偶数 ,记 n = 2m ,k = 0 ,1 , …,m - 1 ,在第 k 个小区间上用三 个节点 (x2k ,f 2k) , (x2k +1 ,f 2k +1) , (x2k +2 ,f 2k +2) 作二次 插值函数 Sk (x) ,然后积分可得 : ∫x2k+1 x2k Sk (x) dx = h 3 [f (x2k) + 4f (x2k +1) + f (x2k +2) ] ,求出 m 段的和就得到整个区间上的近似 积分 Sn ,即 Sn = [f (x0) + f (x2m) + 4 ∑ f (x2k +1) + m- 1 k = 0 h 3 m- 1 2 ∑ k = 1 f (x2k) ] . ④ ④式称为辛普森公式. 3 　求数值积分的 MATLAB 积分命令 3. 1 　矩形公式命令 su m 命令格式 : su m (y) 　% 输出一个向量 y 的分量的和 ,按矩 形公式 ①, ②计算积分的近似值. 3. 2 　梯形公式命令 t rapz 命令格式 : t rapz (x ,y) 　% 输入向量 x = [ x0 ,x1 , …,xn ] , 输出同维数的向量 y = [f0 ,f 1 , …,f n ] ,按梯形公式计 算积分近似值 ; t rapz (y) 　% 按梯形公式计算积分 ,但取步长 h = 1. 3. 3 　辛普森公式命令 quad 命令格式 : quad (′f un′, a , b) 　%计算函数 f un 从 a 到 b 的定积分 ,自动选择步长 ,误差为 10 - 3 ; quad8 (′f un′, a , b) 　% 用高精度计算 ,效果比 quad 更好. 就例 1 中的定积分 ,分别介绍三种命令的用法 , 并对计算结果进行比较. 　　(1) 用左 、右矩形公式命令 ( pi/ 2) / (50 - 1) ( pi/ 2) / (50 - 1) 输入命令 : > > x = linspace (0 , pi/ 2 ,50) ; y = sin ( x) > > I1 = sum ( y (1 :49) ) > > I2 = sum ( y (2 :50) ) 结果为 I1 = 0. 983 9 , I2 = 1. 015 9. (2) 用梯形公式命令 输入命令 : > > I3 = trapz ( x , y) 结果为 I3 = 0. 999 9. 该值就是 I1 和 I2 的平均 值 ,非常接近精确值 1. (3) 用辛普森公式命令 输入命令 : > > I4 = quad (′sin ( x )′,0 , pi/ 2) 结果为 I4 = 1. 000 0. 该值就是精确值 1. 但是 quad 和 quad8 需给出 函数表达式 ,不能像 sum 和 trapz 那样 ,能对用数值 给出的 x , y 进行积分. 4 　数值积分法的应用实例 对于不知道函数具体表达式 , 只知道 x , y 的 一些实测数据的函数 , 可以用数值积分法的梯形公 式求定积分. 以求瑞士国土面积为例说明其应用. 图 1 是瑞士地图 , 为了计算它的国土面积 , 首先对地 图作如下测量 : 以由西向东方向为 x 轴 , 由南向 北方向为 y 轴 , 选择方便的原点 , 并将从最西边 界点到最东边界点在 x 轴上的区间适当地划分为 若干段 , 在每个分点的 y 方向测出南边界点和北 边界点的 y 坐标 y1 和 y2 , 这样就得到表 1 中的数 据 (单位 mm) . 根据地图的比例 ,18 mm 相当于 40 km ,由测量 数据根据下列程序计算瑞士国土的近似面积 : 图 1 　瑞士地图 

06 石家庄职业技术学院学报 第 20 卷 　 表 1 　瑞士地图测量数据 13 47 70 17. 5 50 72 34 50 93 40. 5 44. 5 38 100 30 110 48 30 56 34 110 110 104 106. 5 111. 5 118 123. 5 136. 5 142 7 44 44 96 43 10. 5 45 59 101 37 x y 1 y 2 x y 1 y 2 121 124 121 121 33 28 32 121 65 122 55 116 54 83 52 81 91 46 118 61 36 117 146 50 82 68. 5 76. 5 80. 5 34 118 150 66 86 41 116 157 66 85 45 118 158 68 68 　　x = [ 7 10. 5 13 17. 5 34 40. 5 44. 5 48 56 61 68. 5 76. 5 80. 5 91 96 101 104 106. 5 111. 5 118 123. 5 136. 5 142 146 150 157 158 ] (40/ 18) ; y1 = [ 44 45 47 50 50 38 30 30 34 36 34 41 45 (40/ 46 43 37 33 28 32 65 55 54 52 50 66 66 68 ] 18) ; y2 = [ 44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85 68 ] (40/ 18) ; trapz ( x , y2) - trapz ( x , y1) 运行该程序得结果为 42 414 km2 . 5 　结束语 使用 MA TLAB 命令求积分不仅可对已知解析 表达式的函数进行积分 ,还可以对一些离散点的数 据进行积分 (近似值) ,通过输入简单的命令就可得 到结果 ,轻松完成手工计算难以处理的含有大量数 据的工作. 参考文献 : [ 1 ] 　柯善军. 高等数学与应用实验 [ M ] . 北京 : 北京航空航天大学 [ 2 ] 　萧树铁. 大学数学数学实验 [ M ] . 北京 : 高等教育出版社 , 出版社 ,2007. 1997. 责任编辑 :金 　欣 Methods to solve numerical integration by using MATLAB CHEN Pei ninga , 　L IU Jingb (a. Department of Information Technology ; b. Department of Electric and Electronic Engineering , Shijiazhuang Vocational Technology Institute ,Shijiazhuang ,Hebei 050081 ,China) Abstract :This paper discusses ,wit h illustrations ,t he reasons and met hods of using matlab to solve t he nu merical integration. Key words :MA TLAB ; numerical integration ; rectangular formula ; trapezoid formula ; Simpson formula