数值积分范文

数值积分范文第1篇

关键词： 广义积分 无穷限 无界函数 蒙特卡罗算法

MatLab是一个开放的数学应用软件，由美国的CleveMoler博士研发而成，以矩阵运算为基础，将计算、可视化、程序设计、仿真模拟融合一体，具有工程计算、符号运算、建模仿真、数据分析、图形演示、程序设计等强大功能。MatLab以其强大的功能获得广大科技人员的一致认可，同时也越来越多地被应用在高等数学的教学中，为学生能更好地掌握数学知识、应用数学理论提供了良好的平台。笔者在此试图对MatLab在高等数学数值积分相关内容中的一些使用方法展开讨论，以抛砖引玉，与同仁共勉。

一、对定积分定义的理解

微积分是高等数学中的重要内容，而定积分是其中的重中之重，理解了它对我们以后学习其他相关课程和内容将会带来很大帮助。按照定积分的定义，可将定积分简单解释为被积函数在积分区域被分割成的小积分区间长度与区间上任一点函数值乘积的累加和在积分区间的长度趋于无穷小时的极限，即

2.计算无界函数的广义积分

计算结果如表2。表中对随机点总个数分别为10、100、1000、10000时，做的三次实验结果：

五、小结

在高等数学积分教学中使用MatLab，既提高了学生的学习兴趣，又加深了学生对数学这门枯燥艰深课程内容的理解。同时还能提高教师的素质，也能鼓励学生通过编制一些简单程序来提高动手能力，实在是一件“一箭数雕”的好事。关于这方面的研究，还有待于进一步拓展。

参考文献：

［1］张志涌.精通MATLAB6.5版.北京航空航天大学出版社，2003.3.

［2］胡良剑，孙晓君.MATLAB数学实验.高等教育出版社，2006.6.

［3］陈杰.MATLAB宝典.电子工业出版社，2007.1.

［4］同济大学教研室.高等数学.高等教育出版社，1996.12.

数值积分范文第2篇

【关键词】Newton-Cotes积分公式；matlab；稳定性

1 Newton-Cotes数值积分公式的matlab实现

Newton-Cotes数值积分公式是插值型的，其matlab实现为下面的函数文件：

function y=New_Cotes（a，b，n）。

n=input（'n='）；% n为求积节点的个数

a=input（'a='）；% a为积分下限

b=input（'b='）；% b为积分上限

syms t；

sum=0；

h=（b-a）/n；% h为步长

for i=1：n+1

s=sym（1）；

for j=1：n+1

if j～=i

s=s*（t-j+1）/（i-j）； % 计算连乘

end

end

s（i）=int（s，0，n）； %计算科特斯系数保存在s

y（i）=func（a+（i-1）*h）；%计算函数在节点上的函数值

sum=sum+y（i）*h*s（i）； %计算Newton-Cotes数值积分

end

sum=vpa（sum，6）%显示计算结果，有效数字位数为6。如果用上面的m文件求例1，只需要定义函数func为被积函数，然后运行Newton-Cotes，输入n，a，b，可以得到计算结果。

2 数值算例

例：利用牛顿科特斯公式计算定积分

，取等分区间的份数n=2，4，6，8，10结果如表：

通过上表可以看出当n≤7时，误差能得到有效控制，计算是稳定的；n≥8时，误差不能得到有效控制，计算是不稳定的.同时In（f）也不一定收敛于I（f）。

参考文献：

[1]黄友谦，李岳生.数值逼近（第二版）[M].北京：高等教育出版社，1987

[2]李庆扬，王能超，易大义.数值分析（第三版）[M].武汉：华中科技大学出版社，1986

[3]王能超.数值分析简明教程[M].北京：高等教育出版社，1984

[4]苏金明，阮沈勇.Matlab6.1实用指南[M].北京：电子工业出版社，2002

[5]王晓霞，王治和.求积公式及其误差分析[J].西北师范大学学报，2011，（03）

[6]杨平霞，陈红斌.一类复合插值型求积公式的构造方法[J].宜春学院学报，2011，（8）

资助项目：

四川省高等教育”质量工程”资助。

作者简介：

数值积分范文第3篇

关键词：高斯数值积分方法；积分变量；积分区间；非圆弧拱；多拱梁法

1问题的提出

高斯数值积分方法是一种节点很少、精确度很高的方法，它的特点是节点不等距，计算精度很高，一般利用正交多项式的有关关系式来确定其节点位置和系数.当节点为n时，其代数准确度可达2n-1次.如节点数为3,则求积公式对于任意5(=2×3-1)次多项式都是准确的，这样的精度完全可应用于拱坝程序中的拱段及梁段的计算.笔者曾在圆弧拱多拱梁法程序中，广泛采用3节点高斯数值积分，取得成功[1].但将这一方法推广应用于非圆弧拱多拱梁法程序，却有一定的难度.

由计算数学可知[2]，一般定积分式与高斯积分式的变换形式为：

式中:l=(b-a)/2，为积分限变换系数；n为高斯节点数；ξi为高斯积分节点坐标；gi为与ξi对应的高斯积分系数.

如所周知，拱圈形、载常数计算公式为：

(1)

(2)

式中:S为弧长.

当采用高斯积分公式且为等截面圆拱时，上式相应改为：

(3)

(4)

式中:φ为与弧长S对应的中心角；r为中心半径；E为坝体弹模；Ii,MLi为与高斯节点i对应的截面惯性矩及静定力矩.从上式可见，由于等截面圆拱r为定值，存在简单的dS=rdφ,S=rφ的关系，所以积分变量由S变为φ，积分区间由弧段0～S变为中心角0～φ.在等截面圆拱计算中，许多参数(如坐标及静定力系等)可直接由中心角用显式求得.积分变量由S改为φ后，可大大简化计算.

然而，对于非圆弧拱，问题要复杂得多.因为非圆弧拱的曲率半径和曲率中心处处都在变化，往往不能用简单的显式来表示某一拱段中心角与弧长的关系.并且某些曲线计算弧长也很麻烦，甚至不易用显式求得.如何采用适宜的积分变量和积分区间，是迫切需要解决的问题.

2不同曲线积分变量区间的选取

2.1基本资料5种非圆弧拱示意见图1.几种非圆弧拱曲线方程及有关公式见表1.

图15种非圆弧拱平面示意

表1几种非圆弧拱曲线公式

类型抛物线对数螺旋线椭圆双曲线

曲线方程y=x2/2Rρ=ρ0eaψ

a=cosβ

曲率半径r=(R2+x2)3/2/R2r=Reaψ

r=a2b2x

[(a-y)2/a4+x2/b4]3/2r=a2b2x

[(a+y)2/a4+x2/b4]3/2

弧长S=1/2{xA/cosψA+RLn[(sinψA+1)/cosψA]}S=R/a(eaψA-1)

dSReaψdψ

备注R为拱冠曲率半径β为切线角,R为拱冠曲率半径a为长轴之半,b为短轴之半a为实轴之半,b为虚轴之半

2.2弧的微分公式表1所列弧的微分(dS)算式，除了对数螺旋线稍简单外，其余3种曲线的算式都比较复杂.若从曲线的一般性质来看，当曲线方程y=f(x)时，则弧的微分为：

式中:y′=tgφ，为函数y=f(x)在点x的导数，φ为过点(x,y)的切线与x轴的交角，不难看出，此角与该点的中心角相等.将y′=tgφ代入上式：

所以dS=dx/cosψ(5)

将式(1)作为抛物线、椭圆、双曲线3种曲线弧的微分一般公式，比表1中所列dS算式要简捷得多.式(1)也适用于其他一阶导数存在的任何曲线.

2.3积分变量区间的选取根据各类曲线的性质，选取3种积分变量区间，分述于后.

(1)对于抛物线、椭圆、双曲线3种曲线，采用式(1)所列弧的微分一般公式.此时积分变量为x，积分区间为与S相应的x变化区间.高斯积分时各项均应乘以1/cosφi,φi为与节点i对应的中心角.如A1算式改为:

(6)

(2)对于对数螺旋线，采用表1所列算式：dS=Reaφdφ.此时积分变量为φ，积分区间为与S相应的φ的变化区间，各高斯积分项i均应乘以Reaφi，如A1算式为：

(7)

上式除了高斯积分项乘以eaφi外，其余与圆弧拱相似.

图2五心拱平面示意

(3)对于五心拱，其中弧段为圆弧，算法与等截面圆拱相同，当其边弧段为变截面时，则中心拱弧线为非圆弧曲线，且不便用显式表达.仔细考察该段曲线，发现它与以R3=(RM+RD)/2为半径的圆弧很相近(RM,RD分别为边弧外、内半径)，此圆弧中心O3在边弧起始截面外弧中心O1与内弧中心O2联线中点处，如图2所示.

边弧段计算时，积分变量为φ，积分区间为与边弧S相应的φ的变化区间φ3，如A1算式为：

(8)

上式形式上与圆弧拱一样.

3算例

以上述3种积分变量区间的选取方式，计算各类曲线的半拱弧长，举例于下.

设采用3节点高斯数值积分方法，节点坐标及高斯积分系数列于表2.如以中心角φ为积分变量，以Δφ为积分区间，则与节点i相对应的φi=Δφ(1+ξi)/2.又如以水平坐标x为积分变量，以Δx为积分区间，则与节点i相对应的xi=Δx(1+ξi)/2.

表2节点坐标及高斯积分系数

节点号i节点坐标ξi高斯积分系数gi

1-0.7745966910.555555582

200.888888896

30.745966910.555555582

3.1求抛物线、椭圆、双曲线拱半拱弧长例1：设抛物线拱拱冠曲率半径R=140m，拱端中心角φA=45.32°，拱端坐标xA=141.5726m，求半拱弧长S，可以采用两种方法.一种是按表1所列S的公式直接计算，这是理论积分后的公式，是精确的.另一种是采用高斯数值积分方法，以x为积分变量，xA为积分区间.

(1)按理论公式

S=1/2{xA/cosψA+RLn[(sinψ\-A+1)/coxψA]}

算得S=162.921371m.

(2)按高斯数值积分方法，应有

算得gi/cosφi=2.301610646,S=(141.5726×2.301610646)/2=162.922485m.

该数值积分值与理论计算值162.921371m相比，仅相差0.001114m，相对误差仅为6.8×10-6.

例2：设椭圆拱拱冠曲率半径R=164.51m，拱端中心角φA=51°，坐标xA=141.5843m，长轴之半236.9m，短轴之半197.42m，求半拱弧长.例3：设双曲线拱拱冠曲率半径R=152.71m，拱端中心角φA=41°，坐标xA=141.58421m，实轴之半954.41m，虚轴之半381.77m，求半拱弧长.

上两例因两种曲线无理论积分公式直接用显式计算弧长S，只能与表1中dS算式的高斯数值积分值相比.用dS算式直接数值积分时，积分变量、积分区间与前述方法一样，仍为x及xA，但每一高斯积分项不用除以cosφi，而是乘以dx前的算式等，可见后一算法较繁.两例成果S及比较见表3.

表3椭圆拱、双曲线拱计算成果比较

类别例2椭圆拱例3双曲线拱

原dS算式成果164.143661158.610046

dS=dx/cosψ成果164.143737158.610031

两种算法差值/m0.0000760.000015

相对误差4.63×10-79.457×10-8

由表3可见两种算法成果非常接近.

3.2求对数螺旋线拱半拱弧长例4：设对数螺旋线拱拱冠曲率半径R=154m，拱端中心角φA=46°，切线角β=52°，a=ctgβ=0.781285626,以中心角φ为积分变量，积分区间为拱端中心角φA，求半拱弧长S.

(1)按理论公式计算

S=R/a(eaψA-1)=171.9726625m

(2)按高斯数值积分计算

两种算法成果差值为0.0000061m，相对误差仅为3.547×10-8.

3.3用近似方法求五心拱边例5：设五心拱半拱边弧夹角10°，外半径RM=290m，内半径RD=193.514m，边弧起点拱厚6.466m，拱端厚度8.466m，求半拱边弧长S.

(1)用较精确的计算公式S=φ3×(RA+2×R3)/3.式中3=(RM+RD)/2;φ3为以O3为近似中心的边弧夹角，以弧度计；RA为边弧拱端点与O3联线的长度，见图2.算得RA=241.56728m，φ3=0.206893551弧度，R3=241.757m，从而算出S=50.00488034m.

(2)用近似计算公式S＝φ3×R3=50.01796429m.

两种算法所得的差值为0.013m，相对误差为2.6×10-4.

由以上成果可知，以R3为半径，以φ3为中心角所得的圆弧与边弧的中心弧很近似，因此在高斯数值积分计算时，可以采用较简单的积分变量φ及积分区间φ3.

4结语

综上可知，高斯数值积分应用于非圆弧拱时，需针对非圆弧拱曲线的性质和特点，选取适宜的积分变量与相应的积分区间，这样常可收到事半功倍之效.如对于抛物线拱、椭圆拱、双曲线拱，采用dS=dx/cosφ的通式，既避免了各种曲线的繁复计算，也便于程序规格化.对于对数螺旋线拱，则利用对数函数微分、积分都简单以及极坐标方程弧的微分的特点，积分变量选取φ而不选取x，既大大简化了计算，也很容易利用圆弧拱的算法稍加变换.对于五心拱，在控制误差足够小的前提下，边弧线采用近似圆弧，大大简化了计算.

将3节点高斯数值积分方法推广应用于各种非圆弧拱坝多拱梁法程序，对于提高拱坝程序的计算速度和精度，有很大价值；对于将来推广应用高斯数值积分于各类复杂结构的分析、计算，也有一定的启发和借鉴作用.

参考文献

[1]黎展眉.高斯数值积分方法在多拱梁法程序中的应用[J].砌石坝技术，1986(2).

数值积分范文第4篇

关键词 数值分析 数值积分 算法 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

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

(z>=0 & z

运行结果：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.

数值积分范文第5篇

关键词:数值积分;Richardson外推;奇点;Romberg算法

中图分类号:O241.4

文献标识码:A文章编号:1672-8513(2010)01-0016-04

A Numerical Integration Method for a Class

of the Functions with Singular Points

YANG Lufeng1, JIN Yunchao2

(1.Research Institute of Numerical Computation and Engineering Application, The North University for Ethnics, Yinchuan 750021, China;

2.Guangdong Industry and Trade Vocational School, Foshan 528237, China)

Abstract:

For a class of the functions which contains singular points, the Romberg method is ineffective because of the singularities. This paper adopts the trapezoidal asympotic formula, and provides a novel formulation deduced for the functions, which takes the influence of the singularities into account, like Romberg method. The numerical tests show that the formulation is quite efficient for the class of these functions.

Key words:

numerical integration; Richardson extrapolation; singular points; Romberg method

数值微分与数值积分积分在工程实际中得到广泛应用［1］，国内外众多学者在数值积分领域提出了许多新的计算方法［2-3］,Romberg算法因其收敛速度快被得到广泛应用［4-5］，但对于被积函数存在奇点(被积函数或导数不存在或者不连续点)的一类积分问题，如：I=∫10 xdx和I=∫10 3x•e－xdx等，直接利用Romberg算法求解,计算效果很差. 文中针对I=∫10 xαf(x)dx, 0

1 Richardson外推方法

对于积分问题

I=∫ba f(x)dx,(1)

根据Euler-Maclaurin公式［6］知,对于(1)式所示的积分形式有

∫ba f(x)dx=h2∑2n－1i=0［f(xi)+f(xi+1)］+∑m－1k=1A2kh2k［f(2k－1)(a)－f(2k－1)(b)］－A2m(b－a)h2mf(2m)(ξ).

上式右端第1项恰好为积分I在长度为h=(b－a)2n的子区间上的梯形估计,因此有

I=T(h)+c1h2+c2h4+…+ckh2k+O(h2k+2),(2)

其中ck=f(2k－1)(b)－f(2k－1)(a)为与h无关的常数,满足Romberg公式的外推条件，因此可利用Richardson外推技术得到相应的Romberg算法.但是当被积函数f(x)存在奇点或者其f(2k－1)(b)=f(2k－1)(a)时，利用经典的Romberg算法进行计算时，就会失去或降低外推加速的效果［7］.

定理 若g(x)在［0,1］上充分可微，对于积分I=∫10xαg(x)dx，(0

T(h)=I+c1hβ1+c2hβ2+c3hβ3+….(3)

其中

T(h)=h2［g(0)+2∑n－1i=1g(xi)+g(1)］,xi=ih,h=1n,(β1,β2,β3,β4,β5,…)=(1+α,2,2+α,3+α,4,…).

证明 因为0

令 f(x)=xαg(x),F(x)=∫h0f(t)dt .

按Taylor展开式

f(x0+h)=f(x0)+hf ′(x0)+h22!f ″(x0)+h33!f (x0)+…,

有f(0)=f(h)－hf ′(h)+h22!f ″(h)－h33!f (h)+…,(4)

∫h0f(x)dx=F(h)－F(0)=hF′(h)－h22!F″(h)+h33!F (h)－…=

hf(h)－h22! f ′(h)+h33!f ″(h)－h44!f (h)+….(5)

杨录峰，金云超：一类含奇点函数的数值积分方法

(5)－h2×(4)整理得

∫h0f(x)dx=h2［f(0)+f(h)］－h312f ″(h)+h424f (h)－h580f(4)(h)+….

另一方面，在区间［h,1］被积函数满足Euler－Maclaurin求和公式有

∫1h f(x)dx=h［f(h)+2∑n－1i=2f(xi)+f(1)］－h212［f ′(1)－f ′(h)］+h4720［f (1)－f (h)］－….

两式相加，得

∫10 f(x)dx=T(h)+h212f ′(h)－h312f ″(h)+29h4720f (h)+…－h212f ′(1)－h4720f (h)+….(6)

因为f(x)=xαg(x),0

f(h)=hαg(h)=g(0)hα+g′(0)1!hα+1+g″(0)2!hα+2+… .

所以h2f′(h)=h2{αg(0)hα－1+(α+1)g′(0)1!hα+(α+2)g″(0)2!hα+1+…}=

αg(0)hα+1+(α+1)g′(0)1!hα+2+(α+2)g″(0)2!hα+3+… .

同理 h3f″(h)=α(α－1)g(0)hα+1+(α+1)αg′(0)1!hα+2+(α+2)(α+1)g″(0)2!hα+3+…

带入(6),得

∫10 f(x)dx=T(h)+τ1hα+1+τ2hα+2+…+τ1h2+τ2h4+…

将指数按大小次序排列整理得：

T(h)=∫10 f(x)dx+c1hα+1+c2h2+c3hα+2+c4hα+3+c5h4+…证毕.

将步长h折半带入(3)式,有

T(h2)=I+c1(h2)β1+c2(h2)β2+c3(h2)β3+…,(7)

(7)×2β1－(3),得

2β1T(h2)－T(h)2β1－1=I+γ2hβ2+γ3hβ3+….

通过上面的计算，在没有额外增加计算函数值次数的前提下, 将计算误差从O(h2)提高到O(h4).如此循环下去即可得到此类积分问题的加速算法――类Romberg算法.

2 算法描述

类Romberg算法的计算公式为：

T(0,0)=b－a2［f(a)+f(b)］,

T(n,0)=12T(n－1,0)+b－a2k∑2n－1i=1f(a+(2i－1)b－a2k),T(n,m)=T(n,m－1)+12βm－1［T(n,m－1)－T(n－1,m－1)］. （8）

在实际计算过程中并不能无限制的外推下去，因为,式(8)表明每步外推相当于对当前步的积分近似值加一个修正，当m较大时，修正量很小，几乎起不到修正的作用，并且由于计算量的增加增大了舍入误差的影响.因此，在类Romberg算法的计算中通常设定外推次数为6~8，当超过规定的外推次数，仍没有达到要求的控制精度时，以后的计算只将积分区间进行对分而不再进行进一步外推.

3 数值实验

例1 分别利用Romberg算法和改进类Romberg算法求积分

I1=∫10 x•cos xdx .

计算结果见表1，2种算法对分区间与相邻2次计算的误差关系如图1所示.

其中k表示积分区间对分次数，n=2k表示子区间数，R(n),RI(n)分别表示积分区间分成n个子区间时Romberg算法和改进类Romberg算法计算的积分值,ER,ERI分别表示其相邻2次对分区间积分值的绝对误差.

例2 分别利用Romberg算法和改进类Romberg算法求积分

I2=∫103x•e－xdx.

计算结果见表2，2种算法对分区间与相邻2次计算的误差关系如图2所示

4 结论

对不满足Richardson外推条件的一类含奇点的函数I=∫10xαf(x)dx,0

参考文献:

［1］唐民刚,谢绍龙. 广义CH方程的周期波解及数值模拟［J］.云南民族大学学报:自然科学版,2007, 16(1):13-17.

［2］熊华,杨国孝.一类振荡函数的数值积分方法［J］.北京理工大学学报,1999,19(3):280-284.

［3］ WABG X H,HE Y G,ZENG Z Z.Numerical integration study based on triangle basis neural network algorithm［J］. Journal of Electronics & Information Technology,2004,26(3):394-399.

［4］ALLAHVIRANLOO T.Romberg integration for fuzzy functions［J］.Applied Mathematics and Computation, 2005,168:866-876.

［5］MASTROVIC M,OCVIRK E.An application of Romberg extrapolation on quadrature mathed for solving linear Volterra integral equations of the second kind［J］. Applied Mathematics and Computation,2007,194:389-393.

［6］KINCAID D.数值分析［M］.3版.王国荣,译.北京：机械工业出版社，2005:414-417.

数值积分范文第6篇

【关键词】Flash；数值积分；物理学；课件

0 前言

物理学是一门实验科学，单纯从数学角度来记忆公式是无法深入理解物理规律的。在物理教学过程中，传统的以粉笔和黑板作为媒介的教学方式比较死板，通过借助计算机技术，将物理运动过程制作成多媒体课件，形象、直观的展示出来，能够加深学生对物理规律的感性认识，从而提高教学的效果。

Flash是常用的动画软件之一，具有使用方便、动画效果好的优点，最为重要的是内置一套Action Script编程语言，能够通过程序的方式实现其他软件难以实现的动态效果。在物理学课件的制作过程中，传统静态动画制作方式过程繁琐，而且难以精确的再现物理运动过程。借助Flash提供的Action Script编程语言，制作过程得到了简化，还具有精确、通用性强的优点。

1 Flash中实时展现物理过程的方法

1.1 基于运动方程的方法

Flash采用帧的方式运行，通过改变图形在不同帧的位置、大小等属性来实现动画的效果，属于时间离散的过程，而真实的物理运动过程则是时间连续的。为了在Flash中精确的再现物理运动过程必须对时间连续的过程进行离散化。

物理运动过程可以由运动方程来描述，例如：

匀速直线运动的运动方程为x=x +v t

匀变速直线运动的运动方程为x=x +v t+ a t

抛体运动的运动方程为x=x +v cos（θ ）ty=y +v sin（θ ）t- gt

运动方程直接给出了位置与时间的关系，通过在每帧中使用方程计算出位置坐标就能再现运动过程。

这种方法精准度高，只有计算过程中的舍入误差，且误差不会累积。但该方法必须事先求出运动方程，而且不同场景的运动方程差异极大，所以通用性不是很好。

1.2 采用数值积分的方法

由于基于运动方程的方法不够灵活，通用性差，有必要直接从影响物体运动的物理规律出发，寻找一种通用的方法。

根据牛顿定律可知：物体的运动过程由初始状态（位置、速度）以及受到的力决定，而日常中出现的力可以看成和时间、物置和速度有关的函数，因此可以用以下微分方程来表示物体运动过程。

x″=f（t，x，x′），x（t ）=x ，x′（t ）=v 式1

在物体的初始位置和速度已知的情况下，通过数值积分的方法，计算出下一帧的位置和速度，然后以此类推，也能够再现物体运动过程。这种方法通用性较好，但精度比采用运动方程的方法要差，因为使用数值积分递推计算位置和速度，不仅存在舍入误差，还有数值积分方法带来的截断误差，且误差会累积。不过通过采用高精度的计算方法，误差能做到可以接受的程度。

2 数值积分过程

2.1 欧拉方法

首先将式1改写为以下形式

采用欧拉方法求解上式的过程如下：

2.2 龙格库塔方法

采用龙格库塔方法求解式2的过程如下：

3 实例及性能分析

以斜抛运动为例，其运动过程可由以下微分方程描述：

x″=0，x（t ）=x ，x′（t ）=v cos（θ ）y″=-g，y（t ）=y ，y′（t ）=v sin（θ ）式5

上述式子第一项描述水平方向的运动过程，第二项描述垂直方向的运动过程。

3.1 采用欧拉方法的程序

程序中sx表示水平方向的位置，vx表示水平方向的速度，sy表示垂直方向的位置，vy表示垂直方向的速度，t表示时间，h表示积分步长。方法caculateAccX和caculateAccY用于求取加速度，与式2中的函数v′=f（t，x，v）对应。

3.2 采用龙格库塔方法的程序

以上程序为水平方向的计算过程，垂直方向的计算过程与之类似，程序中的变量和函数与欧拉方法程序的变量和函数相同。可以看到龙格库塔法的计算过程要比欧拉方法复杂，接下来将会对两者的性能进行比对分析。

3.3 性能分析

图1中为取h=0.2s时的运行结果，图中实线为运动方程表示的运动过程，+记号的点序列表示欧拉方法计算结果，×记号的点序列表示龙格库塔方法计算结果。可以看出欧拉方法在初段与运动方程的结果相近，但随着步数增加，误差越来越大，而龙格库塔法的误差几乎可以忽略。

4 结束语

本文介绍了基于数值积分的物理学Flash课件制作方法，给出了采用两种不同数值积分的实现过程，并对两者的性能进行比对，得出结论：欧拉方法计算过程简单，但误差较大，适合在步长较短且运行时间也比较短的场合使用，龙格库塔法计算过程复杂，但误差很小，适合在步长较长且运行时间也比较长的场合使用。

【参考文献】

[1]陈I敏.龙格-库塔法及其Mathematica实现[J].武汉工程职业技术学院学报，2006，18（2）.

[2]黄晓红，胡振华.浅析龙格-库塔方法[J].黑龙江科技信息，2012（23）.

数值积分范文第7篇

关键词：横截面积；矿井巷道；复化梯形公式

DOI：10.16640/ki.37-1222/t.2017.13.210

煤炭工业是国民经济中的基础，它为经济生产提供原料和能源。煤炭工业生产的顺利进行，一定程度上取决于煤炭工业基本建设及开拓延伸工作能否及时、持续地提供生产煤炭的场地。为了将煤从地下采出，需从地表开始，开凿一系列的井筒、硐室及巷道以便到达煤层，这便是矿山基本建设的主体工程。其中涉及到大量的巷道断面的设计，包括巷道尺寸和横截面积的计算。设计的巷道断面直接作为井下巷道施工的依据，也是进行巷道工程概预算的依据。我国煤矿井下使用的巷道断面形状，按其构成的轮廓线可分为圆形类、拱形类、矩形类和梯形类共四大类，在此选择底板水平、两帮垂直、顶板为弧形的半圆拱断面进行数学建模并求其横截面积。

1 问题分析及模型

矿井巷道分为梯形巷道，三心拱巷道，半圆拱巷道等，本案例选择半圆拱巷道模型并求解其面积。

但实际使用这种方法往往有困y，因为大量的被积函数，找不到用初等函数表示的原函数，此时Newton-Leibniz公式不能直接运用，需要用其他有效的数值积分方法求解。在此，运用数值积分方法中的复化梯形求积公式进行求解。

2 复化梯形求积公式原理

在使用牛顿-柯特斯公式时，通过提高阶的途径并不总能取得满意的效果，为了改善求积公式的精度，一种行之有效的方法是复化求积。将积分区间分割为n等份，步长，各节点为。所谓复化求积公式，就是先用低阶的求积公式求得每个子段上的积分值，然后用作为积的近似值。在子区间上使用Newton-Cotes公式，将分割为等份，步长为，节点为记为 ，在上作的阶Newton-Cotes求积公式

3 算法的Matlab实现

3.1 实验数据

半圆拱巷道截面积求法可用数值积分中复化梯形求积公式

求得。用上述的的公式来计算半圆拱巷道截面积，（其中，矩形长，矩形宽）；同时公式（1）中的积分区间，然后将积分区间进行等分，用表示二分次数，即区间等分数，得到个小区间，，每个小区间的长度为，其中。并且按

进行计算。

3.2 Matlab程序代码

function FHQJ

k=2；

m=4；

a=0； % a，b 为区间

b=2；

epsilon=1e-3； % 精确度

fun=@（x.）sqrt（-x.^2+k*x）+m；

n =1；

h=（b-a）/2；

y0=h*（feval（fun，a）+feval（fun，b））；

yiter=y0；

while 1

step =2^（n-1）；

f=sum（feval（fun，a+（1：2：2*step-1）*h））；

y=y0/2+h*f；

if abs（y-y0）

break；

end

h=h/2；

y0=y；

yiter=[yiter，y0]；

n=n+1；

end

yiter

disp（y）；

Error=double（int（'sqrt（-x^2+2*x）+4'，'x'，0，2）-y）； %%与真值误差

disp（Error）；

4 计算结果及分析

4.1 计算结果

本案例设定误差不超过10-3，在Matlab下运行程序后的结果如下（k表示二分次数）：

当运行到k=8， 即时就能满足与真实值误差不超过10-3，此时误差为0.0011。

4.2 结果分析

复化梯形求积公式能够较准确的得到实验结果yiter=9.5696，用在较少的计算量便能够达到预定精度，得到准确值与近似值的绝对误差Error=0.0011，较好的完成了巷道面积的计算问题。

5 总 结

本文建立了矿井巷道截面积的计算模型，根据半圆拱形截面积函数，运用复化梯形求积公式进行计算，并编制基于Matlab的计算程序。根据复化梯形求积公式的原理将所求函数的积分区间分为若干个小的积分区间，先求出每个小积分区间上的积分近似值，然后再将这些近似值加起来就是我们所要求的横截面积的近似值，函数区间所分的小区间的个数越多，计算结果就越精确，其原因就是所分的区间数越多，计算时每个小区所带来的误差就越小，对其求和所带来的总的误差也就越小，所以最后的结果精度就越高。本文算例结果表明该方法具有较高精度、操作简单等优点，且便于程序化。

参考文献：

[1]孙祥，徐流美，吴清.Matlab7.0基础教程[M].北京：清华大学出版社，2005.

[2]薛毅.数值分析与实验[M].北京：北京理工大学出版社，2005.

[3]林柏泉，崔恒信.矿井瓦斯防治理论与技术[M].徐州：中国矿业大学出版社，1998.

[4]煤矿瓦斯治理与应用总体方案[Z].北京：国家安全生产监督管理总局，2005.

[5]汪卉琴，刘目楼.数值分析[M].北京：冶金工业出版社，2004.

[6]李庆扬，王能超，易大义.数值分析[M].武汉：华中科技大学出版社，2006.

[7]肖伟，刘忠，曾新勇等.Matlab程序设计与应用[M].北京：清华大学出版社，2005.

[8]东兆星，吴士良.井巷工程[M].徐州：中国矿业大学出版社，2009.

数值积分范文第8篇

[关键词]定积分 数值积分 常微分方程数值解法

[中图分类号] O172.2 [文献标识码] A [文章编号] 2095-3437（2015）06-0069-02

科学计算被誉为20世纪最重要的科学进步之一。著名的计算物理学家、诺贝尔奖获得者Wilson教授在80年代就指出：“当今科学活动可分为三种：理论，试验和计算”。中国著名的计算数学家石钟慈院士在其2000年的书《第三种科学方法――计算机时代的科学计算》[1]中高度评价了科学计算在现代科技发展与人类社会进步中的重要作用。目前科学计算已经充分融入到各种科学和工程领域当中，造就了一系列新型交叉学科（如计算生物学等）的产生与发展。

然而，目前绝大多数高等院校的数学课程设计，主打课程依旧是《高等数学》等，鲜有院校面对大范围的学生开展类似《计算方法》、《数值分析》的课程教学。然而，对于绝大多数数学学习者来说，他们更希望将数学作为一个工具来解决他们学习与研究中的问题，而且在很多实际问题中，精确解往往很难得到，这让学习数值计算方法具有重要意义。因此，在有限的数学教学中，如何向学生介绍和渗透科学计算的思想尤为重要。

一、基于定积分定义的数值积分

在定积分的计算时，可以通过求出被积函数的原函数求出定积分的精确值。但是在实际应用过程发现，找出一个函数的原函数并非一件易事，许多函数甚至不存在初等函数表示的原函数，例如[2]：

三、结语

在大学数学的教学中，除了要教会学生基本的数学原理之外，还要注重学生使用数学的能力培养。数值计算一方面可以加深学生对数学概念的理解，另一方面，也为他们今后的学习和工作提供一个强有力的工具。

[ 参 考 文 献 ]

[1] 石钟慈.第三种科学方法――计算机时代的科学计算[M].福建：暨南大学出版社，2000.

[2] 同济大学数学系.高等数学（第六版）[M].北京：高等教育出版社，2007.

[3] 赵晶，李宏伟等.工科数学分析[M].武汉：中国地质大学出版社，2010.

[4] 李庆扬，王能超，易大义.数值分析（第五版）[M].北京：清华大学出版社，2008.

[5] A.Quarteroni，R.Sacco，F.Saleri.数值数学（Numerical Mathematics）（影印版）[M].北京：科学出版社，2006.

[收稿时间]2014-12-23

数值积分范文第9篇

关键词：常微分方程初值问题；自适应；数值积分；Matlab

中图分类号：O13文献标识码：A文章编号：1671―1580（2014）02―0148―03

一、引言

对于一阶常微分方程的初值问题

二、基于自适应数值积分的常微分方程数值算法原理

根据上式，可以近似地取Tk+1-Tk作为当前步近似值Tk+1的误差。若预定精度ε满足∫bkakf（x）dx-Tk+1

三、数值算例

用基于三种自适应的积分方法求解x=1时y（1）的数值解结果和误差情况如表1所示，同用等步长h=0.01时的复化梯形公式、4阶的Runge-Kutta方法和4阶Adams显式公式所求的数值解结果进行了比较，如表1所示，相比之下可看出自适应Cotes积分公式和4阶Runge-Kutta方法结果相对最为准确，但前者用时最少。从整体上看基于自适应梯形积分公式的算法所得到的结果明显比复化梯形公式的误差小，运行时间短。相应的基于自适应Simpson积分公式和自适应Cotes积分公式的算法所得到的结果明显比4阶Runge-Kutta方法和4阶Adams显示公式误差要小，运行时间短。同时对比于梯形公式、4阶的Runge-Kutta方法和4阶Adams显式公式三种方法，可以看出计算常微分方程数值解的单步法迭代公式如果收敛阶数越小，程序运行的时间越长，并且误差相对较大。而多步法在相同条件下，4阶Adams显式公式比4阶的Runge-Kutta方法所得数值解误差大。

用基于三种自适应的积分方法求解x=9时y（9）的数值解结果和误差情况如表2所示，同用等步长h=0.01时的复化梯形公式、4阶的Runge-Kutta方法和4阶Adams显式公式所求的数值解结果和误差情况相比，如表2所示，相比之下可看出自适应Cotes积分公式结果也相对最为准确，而且用时最少。从整体上看基于自适应梯形积分公式的算法所得到的结果明显比梯形公式的误差小，运行时间短。相应的基于自适应Simpson积分公式和自适应Cotes积分公式的算法所得到的结果明显比4阶Runge-Kutta方法和4阶Adams显示公式误差要小，运行时间短。同时对比梯形公式、4阶的Runge-Kutta方法和4阶Adams显式公式三种方法，可以看出计算常微分方程数值解的单步法迭代公式如果收敛阶数越小，程序运行的时间越长，误差则相对较大。而多步法在相同条件下，4阶Adams显式公式比4阶的Runge-Kutta方法所得数值解误差大。

基于三种自适应的积分方法选用的节点如图2所示，基于自适应梯形积分公式的算法选用的节点数最多，自适应Cotes积分公式的算法选用的节点数最少。

[参考文献]

[1]李庆扬，王能超，易大义.数值分析（第5版）[M].北京：清华大学出版社，2008.

[2]胡健伟，汤怀民.微分方程数值方法（第2版）[M].北京：科学出版社，2007.

[3]任玉杰.数值分析及其MATLAB实现：MATLAB6.X，7.X版[M].北京：高等教育出版社，2007.

[4]GeraldCF，WheatleyPO.应用数值分析（第7版）[M].白峰杉改编.北京：高等教育出版社，2006.

数值积分范文第10篇

关键词： 神经网络；数值积分；正弦基函数；收敛性

中图分类号：O241.4 文献标识码：A 文章编号：1671－7597（2012）0220009－01

0 引 言

对给定积分求近似值的过程,即用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。常用的数值积分计算方法有Gauss方法、Newton-Cotes方法等。其中Newton-Cotes方法易出现收敛性很难确保，数值很不稳定，计算量较大等问题。

本文提出的数值积分方法，其基本思想是训练正弦函数神经网络权值

值计算。下面讨论算法，并通过实例比较，证明其有效。

1 算法描述

1.1 基于正弦基函数的神经网络模型

假设函数 ，其周期延拓为 ，表示为傅立叶级数：

1.2 正弦基函数的神经网络算法收敛性定理

定理1 设 为学习率， 是隐层神经元的个数，则当 时，神经网络算法收敛。

1.3 正弦基函数的神经网络算法的数值积分定理

1.4 正弦基函数的神经网络算法步骤

2 数值积分实例

例1 文献[3]用梯形法和 方法在积分区间[0,2]上分别计算

算法中，取神经网络隐层神经元个数为 ，学习率 ，训

表1列出了文献[3]和本文算法的结果

解：用 方法计算该积分时遇到了困难[2]，用

3 结论

由实例可知，本文提出算法计算精度相对比其他算法高，并对被积函数要求较低，稳定性高，适应性强，在工程实际中有一定的应用价值。

参考文献：

[1]沈剑华.数值计算基础[M].上海：同济大学出版社，1999:73-109.

[2]王能超.数值分析简明教程[M].北京：高等教育出版社，1997:66-96.

[3]熊华，杨国孝.一类振荡函数的数值积分方法[J].北京理工大学学报，1999，19(3):280-284.

[4]周永权.多项式函数型回归神经网络模型及应用[J]．计机学报,2003，26(9)：1196-1200.

作者简介：