摘要根据《数值分析》课程的实践性特点，在数值积分内容的教学中融入Matlab实现问题的教学，结合问题教学法以提高学生对数值积分方法的理解和应用能力，有利于研究生创新能力的培养。
关键词 数值分析 数值积分 算法 Matlab实现 问题教学法
中图法分类：G643 文献标识码：A
Teaching Research of Matlab Realization of Numerical Integration in
"Numerical Analysis" Curriculum for Postgraduate Students
GE Cishui
(Department of Mathematics & physics, Anhui University of Architecture, Anhui, Hefei 230022)
AbstractBased on the practical characteristics of numerical analysis curriculum, matlab realizations of numerical integrations are introduced in the teaching procedure. Combined with question-based teaching method, it aims to enhance students" understanding and application capability of numerical integration methods, helpful to train their innovation abilities.
Key wordsnumerical analysis; numerical integration; algorithm; matlab realization; question-based teaching method
0 引言
数值分析也称为数值计算方法，是研究用计算机求数学问题近似解的方法、过程及其理论的一个数学分支。数值分析可以说是科学计算的基础和依托，正如我国著名数学家冯康教授所说：数值分析的发展对于提高计算能力的贡献是与新一代计算机的研制同等重要。在我国，几乎所有工科院校硕士研究生都开设了《数值分析》课程。
Matlab是一个功能强大的科学计算平台，它具有强大的数值计算、符号计算和可视化功能，它提供了大量的函数库、工具箱几乎涵盖了所有的工程计算领域，目前Matlab已成为最为普遍的科学计算工具之一。近年来，人们已意识到在《数值分析》课程的课堂教学和实验教学中引入Matlab科学计算软件的重要性，Matlab软件已替代C语言成为辅助数值分析课程教学的首选，事实上，Matlab软件已成为诸多课程，如：自动控制理论、数字信号处理、动态系统仿真等课程的辅助教学工具，另外掌握Matlab软件本身也对我们开展科学研究和解决工程问题至关重要。
数值积分是《数值分析》课程的重要内容之一，但各种研究生用《数值分析》教材讲解数值积分理论与方法时，对数值积分的软件实现则不加介绍或介绍较少，且融入不够，特别对多元函数的数值积分，仅仅以矩形区域上二重积分为例作简单介绍，这远远不能满足工科研究生的学习和应用需要。
本文根据《数值分析》课程的实践性特点，在数值积分的教学中融入Matlab实现问题的教学，结合问题教学法，以提高学生对数值积分方法的理解和应用能力，有利于工科研究生创新能力的培养。
1 数值积分Matlab实现的问题教学法
数值积分计算的问题很多，如振荡积分问题、无界区域上函数积分问题、无界函数积分问题、高维积分问题等等，但从Matlab实现上来说，到目前Matlab系统还没有直接提供一般区域上的三元及三元以上函数的数值积分指令等。
利用问题教学法来讨论和解决数值积分问题的Matlab实现，可以激发学生学习的兴趣，从所求解问题的性质上找原因、从算法上找原因，积极思维，努力给出解决实际问题的方案，提高综合应用来解决实际问题的能力。
下面列举三个数值积分的Matlab实现问题，来说明问题教学法在课堂教学或实验教学中的应用。
1.1 定积分exsin(1000x)dx的计算问题
下面三种方法以及获得的结果
（1）quad(@(x)exp(x).*sin(1000*x),0,pi)
运行该指令后显示：Warning: Maximum function count exceeded; singularity likely.
显示结果： -3.917208719625272
（2）quadl(@(x)exp(x).*sin(1000*x),0,pi)
运行该指令后显示：Warning: Maximum function count exceeded; singularity likely.显示结果：1.722039000823277
（3）quadgk(@(x)exp(x).*sin(1000*x),0,pi)
显示结果：-0.022140670492099
提出问题：为什么上面三种方法获得的结果会不同呢？
问题分析：这是 = 1000的高频振荡积分问题，quad指令和quadl指令采用的算法不大适合求振荡积分，所给结果不正确。由于本定积分的被积函数的原函数是初等函数，所以可用int指令来获得正确结果：
syms x
vpa(int(exp(x)*sin(1000*x),0,pi),16)
显示结果：0.02214067049210878
quadgk有一定的求振荡积分的功能，上述（3）中所给结果精度较高。注意当振荡频率再高时，如计算exsin(1000x)dx，此时若用指令
quadgk(@(x)exp(x).*sin(10000*x),0,pi)，
计算结果也会出现错误，为获得正确的结果，必须设置较高的“MaxIntervalCount”，如采用指令：quadgk(fun,0,pi,"MaxIntervalCount",10000)。由上可知，求高频振荡积分必须选用quadgk指令、设置较高的“MaxIntervalCount”选项值。
1.2 二重积分的计算问题，其中D是由抛物线y = x2，直线x = 10以及x轴所围成的区域
二重积分的计算问题首先要转化为累次积分的计算问题，上述问题可转化为计算，计算此累次积分方法很多，例如下面指令：
（1）dblquad(@(x,y)(x.^2+y).*sin(x+y.^2).*(y>=0 & y<=x.^2), 0, 10, 0, 100, ...
1e-6, @quadl)
显示结果：-70.483695211568843，运行时间：19.962813 秒
（2）quad2d(@(x,y)(x.^2+y).*sin(x+y), 0, 10, 0, @(x)x.^2, "Abstol",1e-6)
显示结果：-70.483662809994698，运行时间：0.485997 秒
（3）quad2dggen(@y,x)(x.^2+y).*sin(x+y),@(x), @(x)x.^2, 0, 10, 1e-6)
显示结果：-70.483662809308356，运行时间：0.105659 秒
问题：为什么上面三种方法所需时间相差数十倍乃至上百倍？
分析：方法（1）是将一般的积分区域“延拓”为矩形区域，利用dblquad指令计算，这样不可避免地进行了很多0乘运算，费时低效；方法（2）是利用Matlab系统quad2d指令，运行效率尚可，它将一般的积分区域映射到矩形区域，然后利用自适应Lobatton算法来计算，有时需要设置较高的“MaxFunEvals”选项值；方法（3）为NIT数值积分工具箱提供的方法，采用了Gauss算法，精度较高，用时最少。
运行时间的差异与指令所采用的算法有关，同时也与所要求的精度有关。若同一指令，要求的精度较高，可通过设置较小的“tol”值或“Abstol”选项值达到，而这时运行时间则会成倍增加。
1.3 三重积分的计算问题，其中由xoy坐标面与旋转抛物面z = 16 - x2 - y2所围成的立体区域
目前Matlab系统尚没有直接提供一般区域上的三元及以上函数的数值积分指令，下面的方法（1）是将一般积分区域“延拓”为长方体区域，然后利用triplequad指令进行计算
（1）triplequad(@(x,y,z)(x.*z+y.^2).*log(1+x.^2+z).*(x.^2+y.^2<=16) .*...
(z>=0 & z<=16-x.^2-y.^2), -4, 4, -4, 4, 0, 16)
运行结果：1.929759987691869e+003，运行时间：234.935393秒
问题：运行时间较长，怎样高效解决一般立体区域上三重积分的计算问题？
分析：方法（1）是将一般积分区域积分“延拓”为长方体区域积分，这样不可避免地引入了很多0乘运算，故费时低效。能否用解决二元函数数值积分问题的思路来解决此三重积分的计算问题呢？效果较方法（1）如何？
为此先将上面三重积分转化为累次积分
然后利用一元函数和二元函数数值积分指令以及Matlab系统提供的函数arrayfun，进行组合来求一般立体区域上三元函数的数值积分，方法如下：
（2）quad2d(@(x,y) arrayfun(@(x,y)quadgk(@(z)(x.*z+y.^2).*log(1+x.^2+z), ...
0,16-x.^2-y.^2),x,y),-4,4,@(x)-sqrt(16-x.^2),@(x)sqrt(16-x.^2))
运行结果：1.930186638624715e+003，运行时间：4.925542秒
（3）quadgk(@(z) arrayfun(@(z)quad2d(@(x,y)(x.*z+y.^2).*log(1+x.^2+z), ... -sqrt(16-z),sqrt(16-z),@(x)-sqrt(16-z-x.^2),@(x)sqrt(16-z-x.^2)),z),0,16)
运行结果：1.930186640541383e+003，运行时间：0.468731秒
方法（2）是利用quadgk+quad2d组合方法，先二重后一重，方法（3）是利用quad2d+quadgk组合方法，先一重后二重，从结果上看这两种方法不但所获得的结果精度较高，而且用时显著减少，其中方法（3）运行效率更高一些。
2 结语
将数值积分的Matlab实现融入《数值分析》课程的教学，通过问题教学法在课堂教学和实验教学中的应用，有助于学生用数值分析的理论和方法分析计算方法所暴露出的问题，找出失败的原因和解决问题的办法，这样既能加深学生对数值积分方法的理解和记忆，又能提高他们解决实际问题的能力，同时还可以节省时间。
工程上涉及数值积分的问题很多，例如：涉及振荡积分、高维积分、无界函数的数值积分等实际工程问题，都可以作为大型作业，布置给学生课后思考解答，充分发挥学生学习的自主性，锻炼学生查阅资料和使用软件帮助文档功能，有利于学生综合解决实际问题的能力和自主创新能力的培养。
基金项目：安徽省高等学校省级教学质量与教学改革工程项目“《数值分析》课程教学内容优化与组合式教学方法的探索研究”[2008jyxm325 ]
参考文献
[1]孙志忠等.数值分析(第2版)[M].南京:东南大学出版社,2006.
[2]颜庆津.数值分析(第3版)[M].北京:北京航空航天大学出版社,2006.
[3]张志涌等.MATLAB R2010a教程[M].北京航空航天大学出版社,2010.
[4]吴晓勤等.数值分析课程中算法设计的教学[J].湖南工业大学学报,2010.24(2):68-71.
“本文中所涉及到的图表、公式、注解等请以PDF格式阅读”
