Matlab中积分函数程序的改进

丁克会;周建来;倪立学

【摘 要】就曲线的长度积分分析了Matlab中三个一维积分函数程序的共性问题,当被积函数可积,函数值在积分区间的端点处无界时,可能会使程序出错,积分失败.分析了程序运行失败的原因,并提出解决方法.提高了Matlab中积分函数的可靠性,增加了用户程序的稳定性. 【期刊名称】《现代机械》 【年(卷),期】2019(000)001 【总页数】3页(P63-65)

【关键词】Matlab;quadl函数;积分;程序 【作 者】丁克会;周建来;倪立学

【作者单位】淮海工学院,江苏连云港222005;淮海工学院,江苏连云港222005;淮海工学院,江苏连云港222005 【正文语种】中 文 【中图分类】TP311

时至今日，Matlab在学校已经成为线性代数，自动控制理论，数字信号处理，动态系统仿真等高级课程的基本教学工具。在设计研究单位和工业部门，Matlab被广泛用于科学研究和解决各种具体问题。它除了具备卓越的数值计算能力外，还有符号计算，文字处理，建模仿真和实时控制等功能。

Matlab的主要特点：1)语言简洁紧凑，库函数极其丰富，使用方便灵活。Matlab编程简单，程序书写形式自由，其库函数都由本领域的专家编写，函数的可靠性高。2)源程序的开放性。除内部函数以外，所有Matlab的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。本文讨论Matlab中几个积分函数的可靠性问题。 1 Matlab中积分函数简介

Matlab的工具箱中有下列几个数值积分函数：quad，quadl，quadv，dblquad，triplequad。其中dblquad和triplequad两个函数通过调用quad，quadv，quadl三个函数来积分。这三个函数要求被积函数连续，当积分区间端点有界时函数很稳定，不会出现问题，当积分区间端点无界而函数可积时，可能求不到积分而出错。因为，当积分区间端点无界时，被积函数的区间端点的函数值为无穷大，这三个函数用同样的方法：通过将计算端点向区间内偏移一微小的距离，重新计算端点函数值，以避开无穷大值(inf)。本文以Matlab 7.1中quadl函数为例来讨论这三个函数的这一共性问题，并约定被积函数为无界函数，在给定区间内函数连续可积。经检验，Matlab 7.1版后续的版本中积分函数存在同样的问题。 2 积分数学模型

本文以椭圆曲线为例，直角坐标表示的标准椭圆曲线方程如下： (1)

其上半部分曲线显式形式为，x∈[-α，α]： (2)

对直角坐标函数曲线y=f(x)，其在[a b]区间的曲线弧长积分表示为：

(3)

对椭圆上半部分曲线，在积分区间内，(3)式显然收敛。

在运行文献[1]中的程序研究圆弧逼近椭圆曲线时，当逼近到曲线的末端时可能出错。出错的可能性与逼近误差的大小和椭圆曲线的参数有关。报错的原因都是积分函数问题。对这一问题，经过单步深入执行程序，发现问题出现在积分计算末端曲线的长上。在逼近计算到曲线末端时，积分区间很短，quadl函数的第73行试图避开区间的末端点函数值为inf(无穷大)时失败，使得整个逼近过程失败。 3 失败原因分析 3.1 机器数系

电子计算机中数的表示大都采用浮点形式表示，机器数系是一个离散的由有限个有理数所组成的集合。如有一台计算机，n位字长(n=3)，采用β进制(β=2)，界码为p，L≤p≤U(硬件常数L=-1，U=2)，则此计算机表示的数的个数是： F=1+2(β-1)βn-1(U-L+1)=33 (4)

该机器数系转换为10进制数为：0，±0.25，±0.3125，±0.375，±0.4375，±0.5，±0.625，±0.75，±0.875，±1，±1.25，±1.5，±1.75，±2，±2.5，±3，±3.5。

可见机器数系是有限不连续，间隔不均匀的。见文献[3]。 3.2 函数eps(n)与eps×n 图1 eps(n)与eps×n函数图像

在Matlab中，eps为浮点相对精度函数。eps为返回从1到下一个最近的较大双精度浮点数的绝对距离，eps=2(-52)。eps(n)是返回从绝对值n到下一个最近的较大浮点数的绝对距离，精度同n。eps(n)与eps×n有着不同的意义，后者表示倍数关系，将它们用图表示，如图1。它们都是非连续的，eps(n)为阶跃平台函数。

当n=2k，k取正整数时，函数值位于某一平台的左端点。如n=16或32时，eps×n=eps(n)。当16≤n<32时，eps(n)=eps(16)，所以，此时有eps(n)3.3 区间端点的取值使端点偏移失败

quadl函数的69和73行，分别表示对函数区间[a b]的左端点和右端点进行回避无穷大处理。73行是对被积函数的右端点进行处理，其语句如下： y(13)=feval(f，b-eps(superiorfloat(a，b))×(b-a)，varargin{：})；

该函数在初始化时，将积分区间分成12段，共13个点，y(13)为区间右端点函数值。当feval(f，b)为inf时，将b用b-eps(superiorfloat(a，b))×(b-a)替换，其中superiorfloat函数返回单精度或双精度字符串，只要a，b有一者为单精度则返回单精度，正常情况下一般都为双精度。而eps(‘double’)=eps，所以b-eps(superiorfloat(a，b))×(b-a)=b-eps×(b-a)。可见函数采用向内偏移eps×(b-a)的距离，来试图避开inf，只有当b-eps×(b-a)≠b时，偏移才成功，才有可能避开inf值。考虑机器数的问题，当b较大，而(b-a)较小时，使得eps(b)与eps×(b-a)不在一个平台上，eps×(b-a)始终在eps(b)平台的左下方，所以b-eps×(b-a)=b，显然，这时不能避开inf值，使得函数运行失败。如23=8是函数eps(n)的一个间断点，9和7在间断点两侧有9-eps(7)=9，增大避开inf的可能性方法可以用b-eps×(b)或b-eps(b)，这两个表达式都有可能偏移成功，而后者较前者偏移的距离一般更小，精度更高，见图1。 3.4 端点偏移成功但函数值仍为inf

上述情况并非对所有情况都适用，就是说：即使区间端点偏移成功，但端点函数值计算后仍然为inf。如椭圆的长短半轴半径分别为40和20，将椭圆的曲线

x2/402+y2/202=1向左平移38，函数区间[-78 2]，即a=-78，b=2，函数如下： (x+38)2/402+y2/202=1

其一阶导数为：

y′=-(x+38)/2/(402-(x+38)2)1/2

y′函数分母含(402-(x+38)2)，当x=b=2，分母为0时，其函数值y′为inf。对自变量x=b-eps×(b-a)，当b=2时有：x=2-eps×(80)，x=2-1.77e-014，可使端点偏移成功。但x值偏移成功后，它与38在进行加运算时丢失精度，变为整数，故偏移失败。

4 函数区间端点初始函数值的改变对积分结果的影响 图2 quadl函数程序框图

quadl函数初始化计算到的13点函数值分11个内点和2个端点，它们的适用范围有所不同，内点的使用仅到判断调整精度结束，见图2第4框。而端点的使用一直持续到积分过程结束。可见，如对初始13点值进行修改，其影响范围是不同的。对内点值的修改仅有可能影响到积分精度，而对端点值的修改肯定会影响到积分结果，两者的影响是两个方面，前者是间接影响，后者是直接影响。若前者有变化，则后者受前者的控制。一般情况下，内点是无需变动的。而端点值的变化是在端点函数值出现inf时才必须改变的。经过实例验证，端点函数值增大时会使得积分值增大，端点函数值减小会使得积分值减小。在积分值是整数时很明显。其积分值的精度是不受影响的。所以调整端点函数值不会改变积分精度，但递归次数会发生变化，数值偏离愈大，函数运行时间会增大。 5 程序改进

考虑以上因素，在函数可积的情况下，当初始化13点函数值的端点值为inf，在偏移失败后，可直接指定端点函数值，以增大函数的可靠性。

以duaql函数为例，在其第73行下插入下两行，用于避开右端点是无穷大的情况。 if ～isfinite(y(13)) y(13)=feval(f，b-eps(b)，varargin{：})；end if ～isfinite(y(13)) y(13)=1e+8；end

在其第69行下插入下两行，用于避开左端点是无穷大的情况。 if ～isfinite(y(1)) y(1)=feval(f，a+eps(a)，varargin{：})；end if ～isfinite(y(1)) y(1)=1e+8；end

程序经过改进后不会再出现上述类似的情况。 6 结语

曲线在端点处的切线垂直，导致其求曲线长的被积函数在端点处无界，当曲线长度较短时，在Matlab中用quadl等函数积分求曲线的长度时可能失败。本文分析了失败原因，提出了改进积分程序的措施，对可积函数在区间端点无界时都有效。提高了积分函数的可靠性，增加了用户程序的稳定性。 参考文献

【相关文献】

[1] 丁克会，席平原，李强化.基于MATLAB单双圆弧混合逼近曲线的算法及实现[J].制造技术与机床，2008(6)：128-131.

[2] 同济大学数学教研室.高等数学[M].北京：高等教育出版社，1996. [3] 孙志忠，袁慰平，闻震初.数值分析[M].南京：东南大学出版社，2002. [4] 李丽，王振领.Matlab工程计算及应用[M].北京：人民邮电出版社，2003.