查看: 15533|回复: 10
请问如何使用EXCEL解决最小二乘估计问题
阅读权限20
在线时间 小时
本帖已被收录到知识树中，索引项：
请问如何使用EXCEL解决最小二乘估计问题
问题如题目,
使用EXCEL解决最小二乘预测计算问题.
我发一个数据.在这个数据基础上加点(X0,Y0)预测..
我经常用这个分析数据.计算很麻烦.在EXCEL上能否完成.
用最小二乘法按拟合曲线，或者用别的非线性回归方法。
最好数学解释和程序都有！感激不尽！
[ 本帖最后由 ccbzj01 于
23:17 编辑 ]
22:56 上传
点击文件名下载附件
3.71 KB, 下载次数: 102
阅读权限20
在线时间 小时
这是数学公式的解决
最小二乘大约是1795年高斯在他那星体运动轨道预报工作中提出的[1]。后来，最小二乘法就成了估计理论的奠基石。由于最小二乘法结构简单，编制程序也不困难，所以它颇受人们重视，应用相当广泛。
如用标准符号，最小二乘估计可被表示为：
& && && && && && && &AX=B& && &(2-43)
上式中的解是最小化 ，通过下式中的伪逆可求得：
& && && && && && &&&A'AX=A'B& &(2-44)
& && && && &(A'A)^(-1)A'AX=(A'A)^(-1)A'B& &(2-45)
& && && && && &&&(A'A)^-1A'A=I (2-46)
& && && && && &&&X=(A'A)^(-1)A'B& & (2-47)
此即最小二乘的一次完成算法，现代的递推算法，更适用于计算机的在线辨识。
最小二乘是一种最基本的辨识方法，但它具有两方面的缺陷[1]：一是当模型噪声是有色噪声时，最小二乘估计不是无偏、一致估计；二是随着数据的增长，将出现所谓的“数据饱和”现象。针对这两个问题，出现了相应的辨识算法，如遗忘因子法、限定记忆法、偏差补偿法、增广最小二乘、广义最小二乘、辅助变量法、二步法及多级最小二乘法等。
阅读权限20
在线时间 小时
& & & & & & & &
这是在某大学论文中发现最好的解决办法.贴上来给大家看看.希望大家能参与讨论.看看有无更好的解决办法.
1、最小二乘法是以误差理论为依据 ,在诸数据处理方法 中 ,误差最小 ,精确性最好。然而在实际教学过程中因其计算比较繁杂 ,学生很少采用这一方法 ,影响了学生运用最小二乘法进行数据处理能力的培养。&&随着计算机的普及,运用最小二乘法进行数据处理有了有力的工具 ,然而采用编写程 序的方法处理数据学生仍感到并不简便。寻找简便易学、容易掌握的计算方法是解决学生掌握最小二乘法进行数据处理的关键。笔者认为运用最常见的学生也比较了解的软件Excel 进行最小二乘法的计算 ,其过程简便而且容易掌握。
2 　运用 Excel 进行最小二乘法的计算
& & Excel&&中有多种工具可用于最小二乘法的计算,其中的“函数”、“图表向导”、“数据分析”在处理数据时各有特点 ,用于最小二乘法计算时不需要编写程序 ,处理数据非常简便。 例:温度变化时 ,测得某铜线的电阻 ,数据记录在 Excel&&中如表1 ,求在0 ℃时铜线的电阻及其温度系数。
& &表 1 　实验数据记录表
& && && && &A& &&&B& &&&C& &&&D& &&&E& && &F& &&&G& &&&H& &&&I& &&&J& &&&K
& & 1& && &t/ ℃&&&&&&4510& &&&&&7010
& & 2& && &R/ Ω&&
　　这一问题可以用 Excel 通过三种不同的方法进行最小二乘法计算。
211 　运用 Excel&&中的“函数”进行计算
& & Excel&&中有各类函数三百余个,分别用于各种不同的计算。其中的线性回归拟合线方
程的斜率函数 SLOPE、线性回归拟合线方程的截距函数 INTERCEPT 以及相关系数函数
CORREL 可用来确定线性方程y = ax + b 的a 、b 两个系数和计算相关系数以判别线性回归
是否合理。
& & 如在上例中 ,在空白的单元格单击“插入”菜单中的“fx 函数”,在弹出的对话框中分别 选中函数名为“SLOPE”、“INTERCEPT”、“CORREL ”的函数,在各自的对话框中输入存放数& &α& &据的单元格区域B2 : K2 和B1 : K1 便可获得斜率 a =& && && &&&R0&&= 0. 00589 ;截距 b =&&R0&&= 1. 433 和& &&&Ω& && && && &α& && &相关系数 R = 0. 9999 的结果。由此可得在0 ℃时铜线的电阻为 11433& && && && && && &&&, 温度系数为& && &=& && &&&-&&3 -&&1 α /&&R = 4. 108&&×10 ℃&&。
212 　运用 Excel&&中的“图表向导”进行
& & “图表向导”是Excel& && &中绘制图表 的工具,提供有十多种“图表类型”。
其中的“XY 散点图”可用来进行回归 分析 ,在生成一张数据分析图时 ,并能方便地得到拟合线方程和相关系数的平方。
& & 单击“插入”菜单中的“图表”,选 中“XY 散点图”;在对话框“步骤之二” 的“数据区域”中输入存放 y&&轴数据的 单元格区域B 2 : K2 ; 在“系列”选项的项
图1 　电阻温度系数的测定& &&&( )
“X&&值 X&&”中输入存放 x&&轴数据的单 元格区域B 1 : K1;在对话框“步骤之三”中确定图形的名称、坐标轴的标题以及网格线, 在 确定图表的插入位置后就完成实验数据分布图。选中所作的图表, 在工具栏单击“图表” 中的“添加趋势线”, 在弹出的对话框选项中“类型”选“线性”; “选项”选中“显示公式”和 “显示R 平方值”的复选框, 便可得到拟合线方程和相关系数的平方。如图 1 所示, 拟合 线方程为 y&&= 0. 0059x&&+ 1. 4333 及相关系数的平方 R2&&= 0. 9999 , 由此也可得在0 ℃时铜线
& && && && & Ω& && && & α& && && && &-&&3 -&&1
的电阻为11433& && &,温度系数&&为4. 108&&×10& &&&℃&&。
213 　运用 Excel&&中的“数据分析”进行计算
& & “数据分析”是Excel&&中为了进行复杂统计或工程分析时节省步骤的一个专用工具。
使用时单击“工具”菜单中的“数据分析”命令。如果“工具”菜单中没有“数据分析”命令 ,& &&&( 则需要安装“分析工具库”。 在“工具”菜单中 ,单击“加载宏”命令 ,在“加载宏”对话框中& &) 选中“分析工具库”。在弹出的“数据分析”对话框中选中“回归”,此工具可通过对一组观 察值使用“最小二乘法”直线拟合 ,进行线形回归分析。在弹出的“回归”对话框“Y 值输入 区域”、“X 值输入区域”中分别输入存放数据的单元格区域 ,选择“输出区域”单选按钮并 输入要显示结果的单元格 ,若选中“线性拟合图”的复选框则可同时生成图表。单击“确 定”就完成了所有计算和作图工作。&&“数据分析”的结果有许多线性回归分析的计算数值。在本例中不但计算出关系数R&&= 0. 9999 ;斜率 a = 0100589;截距 b = 1. 433 , 同时在“标准误差”行中显示测量值yi& && && && && && &的标准 偏差 s (y)& &= 010016 , 在“标准误差”列中显示斜率 a&&的标准偏差s ( a)& && && && && & = 2155&&×10-&&5和截距& && &( )&&b&&的标准偏差s& && &b&&= 0100126 等等分析数据。
214 　各运算方法的特点比较
& &&&从上面 Excel&&中三种计算方法中可看出:利用“函数”运算方法简单 ,但需要记住函数
名称 ,缺点是没有图表显示;利用“图表向导”运算根据对话框 ,所见即所得操作简单 ,数据
和图表都能显示 ,缺点是运算步骤较多;利用“数据分析”运算过程简单 ,运算结果和图表
可一并获得 ,获得的数据分析结果比前两种方法要多而全 ,而过程则简便得多 ,其缺点正
是得到的分析数据太多,许多数据是初步进行最小二乘法计算时所不需要的 ,要能够对太
多的数据进行取舍。对它们各自特点对比如表2 。
& & 表2 　Excel&&中三种最小二乘法计算比较表
& && & 方法& &&&所需步骤& & 获得信息
& && & 函数& &2 步/ 函数& && &截距、斜率、相关系数
& &&&图表向导 共5 步&&图表、拟合线方程、相关系数的平方&&( (&&) (&&) ( ) )& && &数据分析共2 步 截距、斜率、相关系数、标准偏差 s&&y& && &&&、s&&a 、s&&b 、图表、…等
　　根据上面三种计算方法的特点 ,笔者认为在运用最小二乘法处理实验数据时以 Excel 中“数据分析”最为方便又易于掌握,获得的信息也最多。
3 　结束语
& &&&在众多的最小二乘法计算方法中 ,利用 Excel 来处理 ,不用编写程序 ,简便易学 ,容易
掌握,是解决学生掌握最小二乘法进行数据处理的很好工具 ,在培养学生数据处理能力方面很有帮助。
[ 本帖最后由 ccbzj01 于
00:51 编辑 ]
阅读权限20
在线时间 小时
& & & & & & & &
怎么没有人顶啊。
这么好的东西，就是战场上的刺刀和导弹，能解决小问题，也能解决重大。
这用不知道，一用吓一跳！
阅读权限20
在线时间 小时
[em07] ::D :victory:
没人顶,自己顶起来!
EXCEL是个应用工具,在解决实际问题上,要懂得一些专业性的知识,这个问题在科研,质量控制和经济预测方面经常使用.而且准确率很高!
[ 本帖最后由 ccbzj01 于
23:08 编辑 ]
阅读权限20
在线时间 小时
对于相关系数R的分析,可用分析工具.
两个变量间的相关系数可以使用Correl函数,Pearson 和 RSQ函数.
阅读权限90
在线时间 小时
我只想问一下
你附件中的资料
日期和公差存在线性关系嘛？？
阅读权限20
在线时间 小时
需要分析，找出线性关系，我主要做方差分析。回归分析为了预测。
阅读权限10
在线时间 小时
回复 8楼 ccbzj01 的帖子
如果我要拟合不是线，比如一个圆的几个点上，我要拟合圆怎么弄？
阅读权限10
在线时间 小时
研究经济的好工具，最容易使用的偏最小二乘软件PLS+Excel+Word. 网址；
玩命加载中，请稍候
玩命加载中，请稍候
Powered by
本论坛言论纯属发表者个人意见，任何违反国家相关法律的言论，本站将协助国家相关部门追究发言者责任！ & & 本站特聘法律顾问：徐怀玉律师 李志群律师扫扫二维码，随身浏览文档
手机或平板扫扫即可继续访问
力学计算题部分.doc
举报该文档为侵权文档。
举报该文档含有违规或不良信息。
反馈该文档无法正常浏览。
举报该文档为重复文档。
推荐理由：
将文档分享至：
分享完整地址
文档地址：
粘贴到BBS或博客
flash地址：
支持嵌入FLASH地址的网站使用
html代码：
&embed src='/DocinViewer-4.swf' width='100%' height='600' type=application/x-shockwave-flash ALLOWFULLSCREEN='true' ALLOWSCRIPTACCESS='always'&&/embed&
450px*300px480px*400px650px*490px
支持嵌入HTML代码的网站使用
您的内容已经提交成功
您所提交的内容需要审核后才能发布，请您等待！
3秒自动关闭窗口数值积分算法与MATLAB实现+毕业论文60_matlab数值积分-牛bb文章网
数值积分算法与MATLAB实现+毕业论文60_matlab数值积分
话题：，，，
编 号：审定成绩：重庆邮电大学毕业设计（论文）设计（论文）题目： 数值积分算法与MATLAB实现学 院 名 称 ：学 生 姓 名 ： 专 业 ：班 级 ： 学 号 ：指 导 教 师 ： 答辩组 负责人 ：数理学院 数学与应用数学填表时间： 年 月重庆邮电大学教务处制摘 要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。数值积分就是解决此类问题的一种行之有效的方法。积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。【关键词】数值积分 牛顿-科特斯求积公式 高精度求积公式 MATLAB软件ABSTRACTWhen the solution of the definite integral of some function values，because the original function is very complex and difficult to find the elementary function expression, the integral is difficult to accurately calculate, only managed to find the approximate value, and the case is small that allows to direct interface with the Newton - Leibniz formula to calculate the definite integral. Numerical integration is an effective method to solve such problems. The numerical integration is an important branch o therefore, exploring the approximate calculation of the numerical integration method has obvious practical significance. This article departure from the numerical integration problem, described in detail some important numerical integration methods.This paper has introduced detail the Newton - Coates quadrature formula, and in order to improve the calculation accuracy of numerical integration formulas, More precise formulas have Romberg quadrature formulas and the Gauss - Legendre quadrature formula. In addition to the study of these numerical integration algorithm theory, the article also involve what these numerical integration algorithm be programmed by matlab software on the computer, and an example is calculated with a variety of quadrature formulas, finally analysis and comparison to various quadrature formulas calculation error.【Key words】 Numerical integration Newton-Cotes quadrature formula High-precision quadrature formula Matlab software目 录前 言 ........................................................................ 1第一章 牛顿-科特斯求积公式 ................................................... 2第一节 数值求积公式的构造................................................. 2第二节 复化求积公式....................................................... 9第三节 本章小结.......................................................... 12第二章 高精度数值积分算法 ................................................... 13第一节 梯形法的递推...................................................... 13第二节 龙贝格求积公式.................................................... 14第三节 高斯求积公式...................................................... 17第四节 高斯-勒让德求积公式............................................... 19第五节 复化两点高斯-勒让德求积公式....................................... 22第六节 本章小结.......................................................... 23第三章 各种求积公式的MATLAB编程实现与应用 ................................ 24第一节 几个低次牛顿-科特斯求积公式的MATLAB实现 ....................... 24第二节 复化求积公式的MATLAB实现 ...................................... 28第三节 龙贝格求积公式的MATLAB实现 .................................... 33第三节 高斯-勒让德求积公式的MATLAB实现 ............................... 34第五节 各种求积算法的分析比较............................................ 36第六节 本章小结.......................................................... 38结 论 ....................................................................... 39致 谢 ....................................................................... 40参考文献 ..................................................................... 41附 录 ....................................................................... 43一、英文原文 .......................................................... 43二、英文翻译 .......................................................... 52前 言对于定积分?f(x)dx，在求某函数的定积分时，在一定条件下，虽然有牛ab顿-莱布里茨公式I??ba但在很多情f(x)dx?F(b)?F(a)可以计算定积分的值，sinxx况下f(x)的原函数不易求出或非常复杂。被积函数f(x)的原函数很难用初等函数表达出来，例如f(x)?,e?x2等；有的函数f(x)的原函数F(x)存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式。因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值。因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算。而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数f(x)在一些节点上的信息求出定积分的近似值。微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节。数值积分是数学上重要的课题之一，是数值分析中重要的内容之一，也是应用数学研究的重点。随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域。现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义。国内外众多学者在数值积分应用领域也提出了许多新方法。在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等。通过这个课题的研究，我们将会更好地掌握运用数值积分算法求特殊积分函数的定积分的一些基本方法、理论基础；并且通过matlab软件编程的实现，应用于实际生活中。数值积分算法与MATLAB实现+毕业论文60_matlab数值积分第一章 牛顿-科特斯求积公式第一节 数值求积公式的构造大多数实际问题的积分是需要用数值积分方法求出近似结果的。数值积分原则上可以用于计算各种被积函数的定积分，无论被积函数是解析解形式还是数表形式，其基本原理都是用多项式函数近似代替被积函数，用多项式的积分结果近似代替对被积函数的积分。由于所选多项式形式的不同，可以有许多种数值积分方法。而利用插值多项式来构造数值求积公式是最常用的一种方法。 对于积分?ba用一个容易积分的函数?(x)f(x)dx，去代替被积函数f(x)，这样的?(x)自然以多项式Ln(x)为最佳，因为多项式能很好的逼近任何连续函数，而且容易求出其原函数。一、 求积公式的推导在积分区间[a,b]上取有限个点a?x0?x1???xn?b，作f(x)的n次插值n多项式Ln(x)??f(xk)lk(x)，其中，lk(x)(k?0,1,?,n)为n次插值基函数。用Ln(x)近似代替被积函数f(x)，则得?若记Ak?babaf(x)dx??banLn(x)dx??k?0f(xk)?lk(x)dx (1.1)ab?lk(x)dx??ba(x?x0)?(x?xk?1)(x?xk?1)?(x?xn)(xk?x0)?(xk?xk?1)(xk?xk?1)?(xk?xn)(1.2)则得数值求积公式?banf(x)dx??Ak?0kf(xk)(1.3)其中Ak称为求积系数，xk称为求积节点。则称该求积公式为插值型求积公式。知道了插值型求积公式以及其构造方法。为了便于计算与应用，常将积分区间的等分点作为求积节点，这样构造出来的插值型求积公式就称为牛顿-科特斯（Newton-Cotes）求积公式。在积分区间[a,b]上取n?1个等距节点xk?a?kh(k?0,1,?n)，其中h?b?an，做n次拉格朗日插值多项式Ln(x)，因为f(x)?Ln(x)?Rn(x)，所以?baf(x)dx??baLn(x)dx??baRn(x)dx(n?1)?b1?f(x)bl(x)dx??f??k?k??aa??(n?1)!k?0n(?)?n?1(x)dx记Ak??balk(x)dx??ba?n?1(x)(x?xk)?n?1(xk)(1.4)Rn?f??b?(n?1)!n1baf(n?1)(?)?n?1(x)dx (1.5 )截去第二项得?af(x)dx??Ak?0kf(xk)显然Ak与f(x)无关，只与节点xk(k?0,1,?,n)有关。令x?a?th，则当x??a,b?时，t??0,n?，于是?n?1(x)??n?1(a?th)?hn?1t(t?1)(t?2)?(t?n)(1.6)而 ?n?1(xk)?(xk?x0)(xk?x1)?(xk?xk?1)(xk?xk?1)?(xk?xn) ?hnk!(?1)n?k(n?k)! 从而得Ak?记C(n)k?k!(n?k)!(?1)n?kn0(?1)n?khn0t(t?1)??t?(k?1)???t?(k?1)??(t?n)dt?k!(n?k)!n?t(t?1)??t?(k?1)???t?(k?1)??(t?n)dt(1.7)(n)则 Ak?(b?a)Ck故求积公式(1.3)可写成?banf(x)dx?(b?a)?Ckf(xk)k?0(n)(1.8)这就是牛顿-科特斯求积公式，其中Ck称为科特斯系数。 部分科特斯系数取值如下表1.1 科特斯系数具有以下特点[1]n(n)（1） ?Ci(n)?1 （2） Ci(n)?Cn?i(n)（3）当n ? 8 时，出现负数，稳定性得不到保证。而且当n 较大时，由于Runge现象，收敛性也无法保证。故一般不采用高阶的牛顿-科特斯求积公式。 （4）当n ? 7 时，牛顿-科特斯公式是稳定的。表1.1 部分科特斯系数表知道了什么是牛顿-科特斯求积公式，下面我们来看它的误差估计，首先来看看牛顿-科特斯求积公式的截断误差。我们知道牛顿-科特斯求积公式是一个插值型数值求积公式，当用插值多项式Ln(x)代替f(x)进行积分时，其截断误差R?f?即积分真值I和近似值之差，推导如下R?f??I??baf(x)dx??baf(x)dx??baLn(x)dx???f(x)?Labn由插值多项式(x)?dx，的误差估计可知，用n次拉格朗日多项式Ln(x)逼近函数时产生的误差为f(x)?Ln(x)??baf(n?1)(?)(n?1)!n?1(x)dx(1.9)其中?n?1(x)??(x?xi),??(a,b)。对上式两边从a到b作定积分，便可得出0?i?n它的截断误差Rn?f???(n?1)!1baf(n?1)(?)?n?1(x)dx(1.10)二、几个低次牛顿-科特斯求积公式从上面的讨论可知，用多项式近似代替被积函数进行数值积分时，虽然最高次数可是8，但是8次多项式的计算式非常繁杂的。常用的是下面介绍的几种低次多项式。1、 矩形求积公式定义1.1 在牛顿-科特斯求积公式中，如果取n?0，用零次多项式（即常数）代替被积函数，即用矩形面积代替曲边梯形的面积，则有?baf(x)dx??baL0(x)dx?(b?a)c0f(x0)?(b?a)f(x0)(0)(1.11)称式(1.11)为矩形求积公式根据牛顿-科特斯求积公式的误差理论式(1.10)，矩形求积公式的误差估计为R0(f)?f(0?1)(?)(0?1)!?ba?0?1(x)dx?f?(?)(b?a)2、梯形求积公式定义1.2[1] 在牛顿-科特斯求积公式中，如果取n?1，用一次多项式代替被积函数，即用梯形面积代替曲边梯形的面积，则有?baf(x)dx??baL1(x)dx?(b?a)[c0f(x0)?c1f(x1)](1)(1)其中，f(x0)?f(a),f(x1)?f(b)查表可得c0(1)?c1(1)?12代入上式得出?baf(x)dx??baL1(x)dx?(b?a)2)[f(a)?f(b)] (1.12称式(1.12)为梯形求积公式由于用一次多项式L1(x)近似代替被积函数f(x)，所以它的精度是1。也就是说，只有当被积函数是一次多项式时，梯形求积公式才是准确的。根据牛顿-科特斯求积公式的误差理论式(1.10)，梯形求积公式的误差估计为R1(f)?f(2)(?)(2)!?ba(x?a)(x?b)dx??(b?a)123f??(?)f??(?)是被积函数f(x)二阶导数在x??点的取值，??[a,b]3、辛浦生求积公式定义1.3ba[2]在牛顿-科特斯求积公式中，如果取n?2，用二次多项式代替ba被积函数，即曲边用抛物线代替，则有?上式得出f(x)dx??L2(x)dx?(b?a)[c0f(x0)?c1f(x1)?c2f(x2)]a?b2(2)(2)(2)其中，x0?a，x1?baba，x2?b,查表可得c0(2)?c2(2)?1616，c1(2)?23，代入?f(x)dx??L2(x)dx?(b?a)[16f(a)?23f(a?b2)?f(b)] (1.13 )称式(1.13)为辛浦生求积公式，也称抛物线求积公式。它的几何意义是：用过3个点(a,f(a))，(a?b2,f(a?b2))，(b,f(b))的抛物线和x?a，x?b构成的曲边梯形面积，近似地代替了被积函数f(x)形成的曲边和x?a，x?b构成的曲线梯形面积。下面对辛浦生求积公式的误差进行估计。由于辛浦生求积公式是用二次多项式逼近被积函数推得的，原则上它的代数精度为2.但因多项式次数是偶数，数值积分算法与MATLAB实现+毕业论文60_matlab数值积分根据定理1.1可知，它的代数精度为3过(a,f(a))，(a?b2,f(a?b2))和(b,f(b))3个点，构造一个f(x)的三次a?b2)?f?(a?b2)Lagrange插值多项式L3(x),且使L3?(定理得f(x)?L3(x)?f(4)。根据Lagrange插值余项(?)4!?4(x)?f(4)(?)4!(x?a)(x?a?b2)(x?b) ??[a,b]2对上式两边从a到b进行积分，即可得到R2(f)??ba[f(x)?L3(x)]dx?14!?baf(4)(?)(x?a)(x?a?b2)(x?b)dx (1.14)2根据定积分中值定理可知，在[a,b]上总有一点满足下述关系：?baf(4)(?)(x?a)(x?a?b2?f(4))(x?b)dx (?)?(x?a)(x?ab2a?b2)(x?b)dx2通过变量代换t?x?a,dt?dx，很容易求得?ba(x?a)(x?a?b2)(x?b)dx??2(b?a)1205把这个结果代入式(1.14)，便得出辛浦生求积公式的误差估计式R2(f)?14![?(b?a)1205f(4)(?)]??(b?a)28805f(4)(1.15 )(?) ??[a,b]4、科特斯求积公式由于n?3和n?2时具有相同的迭代精度，但是n?2时计算量小，故n?3的Newton-Cotes积分公式用的很少。定义 1.4 在牛顿-科特斯求积公式中，如果取n?4时，牛顿―科特斯公式为?baf(x)dx?b?a90)?7f(x0)?32f(x1)?12f(x2)?32f(x3)?7f(x4)? (1.16称式(1.16)为科特斯求积公式。同理根据式(1.10)可求得其误差估计式R4[f]??2(b?a)b?a6()f9454(6)(?) (??(a,b)) (1.17)三、求积公式的代数精度如果被积函数f(x)为任意一个次数不高于n次的多项式时，数值求积公式一般形式的截断误差R?f??0；而当它是(n?1)次多项式时，R?f??0，则说明数值求积公式具有n次代数精度。一个数值求积公式的代数精度越高，表示用它求数值积分时所需逼近被积函数的多项式次数越高。定理 1.1[3] 牛顿-科特斯求积公式的代数精度等于n，当n为偶数时，牛顿-科特斯求积公式的代数精度等于n?1证明 如果被积函数f(x)是一个不大于n次的多项式，则f(n?1)(x)?0，即R?f??0；而当时任意一个(n?1)次多项式时，f(n?1)(x)?0，故R?f??0 所以，按照代数精度的定义可知，一般情况下，牛顿-科特斯求积公式的代数精度等于n。当f(x)为n次多项式时，f(n?1)(?)?0 ??(a,b)牛顿-科特斯求积公式的代数精度至少等于n。若设是一个(n?1)次多项式，这时f(n?1)(?)为一常数，而R?f???ba(f(x)?Ln(x))dx?bf(n?1)(?)(n?1)!?ba(x?x0)(x?x1)?(x?xn)dx因此，只要证明在n是偶数时，?(x?x0)(x?x1)?(x?xn)dx?0，R?f??0，a定理就可得证。为此，设xi?1?xi?h,(i?0,1,2?n)，令x?a?th,t??0,n?,dx?hdt，于是?(x?x0)(x?x1)?(x?xn)dx?habn?1?nt(t?1)(t?2)?(t?n)hdt由于n为偶数，不妨设n?2k，k为正整数，则t??0,2k?，于是?nt(t?1)(t?2)?(t?n)hdt??2kt(t?1)?(t?k)(t?k?1)?(t?2k?1)(t?2k)dt再引进变换u?t?k，则t?u?k，du?dt，u???k,k?，代入上式右侧，得出?nt(t?1)(t?2)?(t?n)dtk?k??(u?k)(u?k?1)?(u?1)u(u?1)?(u?k?1)(u?k)du??k?ku(u?1)?(u?(k?1))(u?k)du22222最后的积分中被积函数是奇函数，所以积分结果等于零，定理1.1得证。第二节 复化求积公式前面导出的误差估计式表明，用牛顿-科特斯公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点，亦即提高插值多项式的次数。龙格现象表明，这样做并不一定能提高精度。理论上已经证明，当n??时，牛顿-科特斯公式所求得的近似值不一定收敛于积分的准确值，而且随着n的增大，牛顿-科特斯公式是不稳定的。因此，实际中不常用高阶牛顿-科特斯公式。为了提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。一、复化梯形求积公式在实际应用中，若将积分区间分成若干个小区间，在各个小区间上采用低次的求积分式（梯形公式或辛浦生公式），然后再利用积分的区间可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化积分公式的基本思想。以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间的长度趋于零时，小梯形面积之和即就趋于曲边梯形面积的准确值，即定积分的准确值。定义 1.5[4] 将积分区间[a,b]进行N等分，记为h?b?aN,xk?a?kh在每个小区间[xk,xk?1] (k?0,1,?,N?1)上用梯形公式求和，得?baN?1f(x)dx???k?0xk?1xkN?1f(x)dx??k?0h2[f(xk)?f(xk?1)]若将所得的近似值记为TN，整理得?baf(x)dx?h2N?1[fa(?)fb(?)?2fxk(?T)N] (1.18)k?1称式(1.18)为复化梯形公式。记为 TN当N??时， TN?1212N?1N[?f(xk)h?k?0b?k?1baf(xk)h]?ba即TN收敛于?f(x)dxab[?f(x)dx?a?f(x)dx]??f(x)dx如果f(x)?C(2)[a,b]，则在小区间[xk,xk?1]上，梯形公式的截断误差为?xk?1xkf(x)dx?h2[f(xk)?f(xk?1)]??h312f??(?k) ?k?(xk,x?k1)因此RT(f)??baf(x)dx?TN???12h3N?1f??(?k)k?0因为f??(?)区间[a,b]上连续，由介值定理知存在?k?[a,b]，使得f??(?)?1NN?1?k?0f??(?k)从而有RT(f)??baf(x)dx?TN??h312Nf??(?)??(b?a)122)hf??(?) (1.19这就是复化梯形公式的截断误差。下面讨论复化梯形公式的数值稳定性。设计算函数值f(xk)时产生的误差为?k(k?0,1,?,N)，则用式(1.18)计算结果的误差为??h2?0??N?2??k?nhmax??kk?10?k?NN?1?(b?a)max???0?k?Nk因此，无论N为多大，复化梯形公式都是稳定的。二、复化辛浦生求积公式如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就可导出复化Simpson公式。定义1.6[5] 将积分区间[a,b]分成N?2m等分，分点为xk?a?kh，b?a(k?0,1,?,N?1) h?N在每个小区间[x2k?2,x2k]k(?1,2?,N,?上1)。用Simpson公式求积分，则有?x2kx2k?2f(x)dx?x2k?x2k?26h3[f(x2k?2[f(x2k?2)?4f(x2k?1)?f(x2k)])?4f(x2k??)f(xk2?1)]求和得?bamf(x)dx???k?1mx2kx2k?2f(x)dx??3?f(xk?1bah2k?2)?4f(x2k?1)?f(x2k)?整理后得到SN??f(x)dx?h3m?1m[f(a)?f(b)?2?f(x2k)?4?f(x2k?1)] (1.20 )k?1k?1式(1.20)就称为复化辛浦生求积公式。记为 SN如果f(x)?C(4)[a,b]， 则由Simpson插值余项公式可得复化公式的截断误差为RS(f)??baf(x)dx?(2h)5h3m?1m[f(a)?f(b)?2?f(x2k)?4?f(x2k?1)]k?1k?1m???2880k?1f(4)(?) ??[x2k?2,xk2]因为ff(4)(4)(?)为连续，故存在??[a,b], 使得(?)??m1mf(4)(?k)k?1数值积分算法与MATLAB实现+毕业论文60_matlab数值积分代入上式得mRS(f)??k?1?(2h)52880(4)(?)??b?a180hf4(4)(?) ??(a,b) (1.21)式(1.21)表明，步长越小，截断误差越小。与复化梯形公式的分析相类似，可以证明，当N?2m??时，用复化Simpson公式所求得的近似值收敛于积分值，而且算法具有数值稳定性。三、复化科特斯求积公式定义1.7 将积分区间[a,b]等分为N个子区间[x4k?4,x4k]，每个子区间的中点x4k?2,(k?0,1,?,N?1)，子区间长度h?斯公式求和，得b?aN， 在每个子区间上用科特?baf(x)dx?h90NN[7f(a)?32?f(x4k?3)?12?f(x4k?2)?k?1N?1k?1(1.22)h4N32?f(x4k?1)?14?f(x4k)?7f(b)]?CNk?1k?1b?aN式(1.22)就称为复化科特斯求积公式，记为CN ，式中，h? ，xk?a?k(k?0,1,2,?,2N?1)类似地可以推出复化科特斯公式的截断误差为R4(N)(f)??2(b?a)h6()f9454(6)(?) (??(a,b))3) (1.2第三节 本章小结本章节开篇介绍了数值求积公式的构造，主要是用运用插值多项式。接着介绍了几个低次的牛顿-科特斯求积公式，即梯形公式、辛浦生公式、科特斯公式，以及牛顿-科特斯求积公式的改进复化求积公式，并对各个求积公式进行了相应的误差分析。第二章 高精度数值积分算法复化求积公式是提高精确度的一种行之有效的方法，但是在使用复化型求积公式之前必须先给出步长。步长太大精度难以保证，步长太小则又会导致计算量的增加，而事先给出一个合适的步长往往是困难的。在实际计算中常常采用变步长的计算方法，即在步长逐次减半的过程中，反复利用复化求积公式进行计算，并同时查看相继两次计算结果的误差是否达到要求，直到所求得的积分值满足要求为止。下面以梯形公式为例第一节 梯形法的递推在变步长的过程中探讨梯形法的计算规律如下：设将积分区间?a,b?分为N等分，则一共有N?1个等分点xk?a?kh,h?b?aN，k?0,1?N，这里用TN表示复化梯形法求得的积分值，其下标N表示等分数。由余项公式(1.19)可知，积分值为I?TN?b?ab?a2()f??(?1) (a??1?b) 12N再将各子区间分半，使得区间成2N等分。此时所得积分近似值记为T2N，则再由余项公式(1.19)可知，积分值为I?T2N?b?ab?a2()f??(?2) (a??2?b) 122N假定f??(x)在?a,b?上变化不大，即有f??(?1)?f??(?2)，于是得I?TNI?T2N?4，左式也可以写成为I?T2N?13(T2N?TN)?T2N?14?1(T2N?TN)1(2.1)这说明用T2N作为积分I的近似值时，其误差近似为(T2N?TN)。计算过程中3常用(T2N?TN)??是否满足作为控制计算精度的条件。如果满足，则取T2N作为I的近似值；如果不满足，则再将区间分半，直到满足要求为止。 实际计算中的递推公式为T1?b?a2?f(a)?12TN?f(b)?T2N?b?a2NN?j?1f(a?(2j?1)b?a2N)(N?2k?1,k?1,2?)(2.2 )在给定控制参数?后，当满足(T2N?TN)??时，则以T2N作为积分I的近似值。 通过类似的推导，还可以得到下面的结论对于辛浦生公式，假定f(4)(x)在?a,b?上变化不大，则有I?S2N?115(S2N?SN)?S2N?14?12(S2N?SN)(2.3)对于科特斯公式，假定f(6)(x)在?a,b?上变化不大，则有I?C2N?163(C2N?CN)?C2N?14?13(C2N?CN)(2.4 )第二节 龙贝格求积公式梯形法的算法简单，单精度低，收敛的速度缓慢。如何提高收敛速度以节省计算量，这是人们极为关心的课题。由此引出了龙贝格公式。由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果T2N作为积分I的近似值，其误差近似值为(T2N?TN)。可以设想，如果用这个误差作为T2N的31一种补偿，即将T2N?13(T2N?TN)=4T2N?TN4?1作为积分的近似值，可望提高其精确度。直接根据复化求积公式，不难验证SN?T2N?13(T2N?TN)?4T2N?TN4?1(2.5)这说明，将区间对分前后两次复化梯形公式的值，按式(2.1)作线性组合恰好等于复化辛浦生公式的值SN，它比T2N更接近于近似值。同样，根据式(2.3)用S2N于SN作线性组合会得到比S2N更精确的值，且通过直接验证可得CN?S2N?115(S2N?SN)?4S2N?SN4?122(2.6)再由式(2.4)用C2N与CN作线性组合，又可得到比C2N更精确的值，通常记为RN,即RN?C2N?163(C2N?CN)?4C2N?CN4?133(2.7)式(2.7)就称为龙贝格求积公式。上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推方法。我们将序列?TN?,?SN??CN?和?RN?分别称为梯形序列、辛浦生序列、科特斯序列和龙贝格序列。由龙贝格序列当然还可以继续进行外推，得到新的求积序列。但由于在新的求积序列中，其线性组合的系数分别为14m[2]4mm4?1?1与?1?0(m?4)。因此，新的求积序列与前一个序列结果相差不大。故通常外推到龙贝格序列为止。利用龙贝格序列求积的算法称为龙贝格算法。这种算法具有占用内存少、精确度高的优点。因此，成为实际中常用的求积方法。 下面给出龙贝格求积算法的计算步骤：第1步：算出f(a)和f(b)的值，根据公式(2.2)计算T1； 第2步：将?a,b?分半，算出f(公式(2.5)计算S1；第3步：再将区间分半，算出f(a?b?a4)及f(a?3?b?a4)的值，并根据公a?b2)后，根据公式(2.2) 计算T2,再根据式(2.2)和(2.5)计算出T4及S2，再由公式(2.6)计算出C1；第4步：将区间再次分半，计算T8，S4，C2，并由公式(2.7) 计算R1； 第5步：将区间再次分半，类似上述过程计算T16，S8，C4，R2。 重复上述过程可计算得到R1，R2，R4?，一直算到龙贝格序列中前后两项之差的绝对值不超过给定的误差限为止。定义2.1[6] 上述用若干个积分近似值算出更精确的积分近似值的方法，称之为外推法。上述计算步骤也可用表 2.1表示表 2.1 龙贝格计算步骤表数值积分算法与MATLAB实现+毕业论文60_matlab数值积分第三节 高斯求积公式前面介绍的n?1个节点的 Newton -Cotes求积公式，其特征是节点是等距的。这种特点使得求积公式便于构造，复化求积公式易于形成。但同时也限制了公式的精度。n是偶数时，代数精度为n?1，n是奇数时，代数精度为n；我们知道n?1个节点的插值型求积公式的代数精确度不低于n。设想：能不能在区间?a,b?上适当选择个节点x0,x1,x2?xn 使插值求积公式的代数精度高于n？答案是肯定的，适当选择节点，可使公式的精度最高达到2n?1，这就是本节所要介绍的高斯型求积公式。 例：?f(x)dx?A0f(x0)?A1f(x1)?11其中，x0,x1固定在?1,1，A0,A1可以适当选取，只有两个自由度，得到的是梯形公式，其代数精度只有1。如果对求积节点x0,x1也进行适当的选取，将有四个自由度，得到如下公式：?1?1f(x)dx?f(??f这个积分公式的代数精确度为3，这就是高斯型求积公式，上面的求积节点?称为高斯点。一、高斯型求积公式和高斯点对于含有2n?2个参数，AK,xk(k?0,1,?,n)的求积公式：?banf(x)dx??Ak?0Kf(xk)，适当选取这2n?2个参数，可以使得数值积分公式的代数精确度达到2n?1，我们称这一类求积公式为高斯型求积公式，称这类求积公式的积分节点为高斯点，系数AK称为高斯系数。定义 2.2[7] 如果n?1个求积节点的求积公式的代数精度为2n?1，则这n?1个求积节点称为高斯点。因为高斯求积公式也是插值型求积公式，故有结论 ：n?1个节点的插值型求积公式的代数精度d，满足n?d?2n?1二、高斯点的特征定理 2.1[7] 设x0,x1,?,xn是n?1个相异点，以这n?1个点为零点的n?1次多项式为?n?1(x)?(x?x0)(x?x1)则x0,x1,?,xn是高斯点的充要条?(x?xn)件是对于任意不超过n次的多项式q(x)，成立?q(x)?n?1(x)dx?0ab证明 1）必要性。设x0,x1,?,xn为高斯点，则对任意不超过次2n?1的多项式f(x)，均有?f(x)dx??Akf(xk)，则对于任意次数不超过n次多项式q(x)，ak?0bnq(x)?n?1(x)是次数不超过2n?1的多项式，且注意到?n?1(xk)?0，(k?0,1,?,n)，故?banq(x)?n?1(x)dx??Aq(xkk?0k)?n?1(xk)?02）充分性。设对于一切次数不超过n次的多项式q(x)，成立?baq(x)?n?1(x)dx?0，又设f(x)是次数不超过2n?1次的多项式，用?n?1(x)去除f(x)，商q(x)，余r(x)，即f(x)?q(x)?n?1(x)?r(x)，可知q(x)和r(x)均是不超过n次的多项式，从而?f(x)dx?ab?baq(x)?n?1(x)dx??bar(x)dx??bar(x)dx，又因求积公式是插值多项式的构造导出的，由A0,A1,?,An的选取，其代数精确度可以达到n，而r(x)是次数不超过n次的多项式，因此成立?bar(x)dx??Akr(xk)，因?n?1(xk)?0，所以f(xk)?r(xk)，k?0nk?0,1,?,n，故而?banr(x)dx??Ar(xkk?0k)也即?f(x)dx?abn?Ar(xkk?0k)由于f(x)是次数不超过2n?1次的多项式，因此该积分公式的代数精确度至少为2n?1，因而由定义2.2知节点x0,x1,?,xn是高斯点。定义2.3 对于关系?q(x)?n?1(x)dx?0，我们称之为正交性，即?n?1(x)与ab任意n次多项式正交，而这样的多项式类称为正交多项式。三、高斯型求积公式的余项Rn??banf(x)dx??Akf(xk)?k?01(2n?2)!f(2n?2)(?)??ab2n?1(x)dx其中，??(a,b)，?n?1(x)?(x?x0)(x?x1)?(x?xn)。积分，可以计算得到Rn?22n?3[(n?1)!](b?a)34[(2n?2)!](2n?3)f(2n?2)(?) (2.8)第四节 高斯-勒让德求积公式常用的高斯型求积公式有高斯-勒让德公式、高斯-切比雪夫公式、高斯-拉盖尔公式、高斯-埃尔米特公式，下面着重介绍高斯-勒让德公式一、Legendre多项式Pn(x)?1ndnn2?2!dx[(x?1)]，n?0,1,2,? (2.9 )2n称式(2.9)为勒让德（Legendre）多项式，其具有前面提到的正交性性质，即对于任意次数不超过n的多项式q(x)，成立?q(x)pn?1(x)dx?0?11因此，多项式pn?1(x)的零点就是相应的高斯-勒让德求积公式的高斯点。勒让德多项式的前几项如下：p0(x)?1，p1(x)?1d2dx(x?1)?x，p2(x)?212d222?2!dx2[(x?1)]?2212(3x?1)2p3(x)?1218(5x?3x)，p4(x)?318(35x?30x?3)4p5(x)?(65x?70x?15x),?53勒让德多项式的性质：性质1 勒让德多项式的首项系数为 an?(2n)!2(n!)n2性质 2 当n为奇数时，pn(x)为奇函数；当n为偶数，pn(x)为偶函数x)p(x)dx性质 3 对一切次数不高于n?1次的多项式q(x)，有?q(n?110?二、高斯-勒让德求积公式1）当n?0时，p1(x)?x，其零点为x0?0，易得A0?2，?f(x)dx的高?11斯-勒让德求积公式是?f(x)dx?2f(0)，其代数精确度为1?112）当n?1时，p2(x)?112(3x?1)，其零点为x0??2，x1?，此时A0?1,A1?1，?f(x)dx的高斯-勒让德求积公式是：?1?1?1f(x)dx?f(??f，其代数精确度为33）当n?2时，p3(x)?12(5x?3x)，其零点为30，设高斯-勒让德求积公式为：?1?1f(x)dx?A0f(?A1f(0)?A2f令f(x)?1，得到方程A0?A1?A2?2 ○1 令f(x)?x，得到方程?A0?A2?0 ○2 令f(x)?x2，得到方程(A0?A2)?5323○359联立方程○1、○2、○3可解得A0?A2?所以?f(x)dx??11，A1?891??5f(9???8f(0)?5f??，其代数精确度为54）当n?3时，p4(x)?3518(35x?30x?3)，由p4(x)?0，得42x?2，故四个零点为x?，即x0?，x1?x2?，x3?相应的，可以解出A0?A3?12?，A1?A2?12?，高斯-勒让德求积公式为?1?1f(x)dx?A0f(x0)?A1f(x1)?A2f(x2)?A3f(x3)其代数精确度为7综上，高斯-勒让德求积公式的积分节点和积分系数表如 表2.2表2.2 积分节点和积分系数表数值积分算法与MATLAB实现+毕业论文60_matlab数值积分三、高斯-勒让德求积公式的余项对n阶Gauss?Legendre求积公式（共有n?1个求积节点）Rn?1(2n?2)!f(2n?2)(?)???112n?1(x)dx (2.10)其中，??(?1,1)，?n?1(x)?(x?x0)(x?x1)?(x?xn)x0,x1,?,xn为多项式pn?1(x)的零点。积分，可以计算得到： Rn?22n?3[(n?1)!]34[(2n?2)!](2n?3)f(2n?2)(?) (2.11 )四、一般区间上的高斯-勒让德求积公式前面讲解的高斯-勒让德求积公式是计算标准区间[?1,1]上的定积分，对于一般的有界区间[a,b]上的定积分，可以通过变量代换转化为区间[?1,1]的定积分。即在积分?f(x)dx中，令x?abb?a2t?a?b2，当x?a时，t??1；当x?b时，?2b?ax?b?ab?at?1；且dx?bb?a2dt。这个变换为线性变换，其反变换为tb?a2于是成立?af(x)dx??1?1f(b?a2t?a?b2)dt然后，我们可以用[?1,1]区间上的高斯-勒让德求积公式。第五节 复化两点高斯-勒让德求积公式本节对两点Gauss?Legendre公式?1?1f(x)dx?f(?3?f3(2.12)进行推广和复化，得到了新的数值积分公式?1?1f(x)dx?hn?1?2i?0ba(f(xi?12?6)?f(xi?12?6)) (2.13)同时分析了它的积分误差和收敛阶定义 2.4 设I?[8]?若当h?0时成立 f(x)dx，In是一种复化求积公式，I?Inhhlimh?0?C,C?0 (2.14 )则称求积公式In是p阶收敛的。例如，复化的Trapezoid公式和复化的Simpson公式分别具有二阶和四阶收敛性。定理 2.2[8] 设f(x)?C(4)[a,b]，h?公式为b?an，则复化两点高斯-勒让德求积?baf(x)dx?hn?1?24(f(xi?12?6)?f(xi?12?6)) (2.15 )k?0的余项表达式为R[f]?(b?a)h52?135f(4)) (?),??[a,b] (2.16该方法具有四阶收敛性。考虑积分I??1e?x2dx，准确解为 0.?用复化两点高斯-勒让德公式可求得I=0. ，其绝对误差为R=0.，与复化梯形公式，复化Simpson公式以及原两点Gauss-Legendre公式相比，这里构造的复化两点高斯-勒让德求积公式方法在运算量和精度方面都有显著的改善。第六节 本章小结本章主要介绍了精确度比较高的两种数值求积公式，即龙贝格求积公式和高斯型求积公式，对其进行了相应的误差分析。其中高斯型求积公式主要介绍了高斯-勒让德求积公式，并对两点高斯-勒让德求积公式进行了复化。第三章 各种求积公式的MATLAB编程实现与应用MATLAB是由MathWorks公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过MATLAB软件编程实现，以下程序均用MATLAB7.0编写，运行坏境：1、硬件环境CPU（intel Core i3-GHz），内存（2GB昱联）,2、软件环境windows7（32位）操作系统。以下总共编写了六个算法程序，部分代码参考文献[10-13],为了体现程序的正确性，以下程序都以I??10sinxx为例进行运算。原积分的精确值为dx?0.183I??10sinxx第一节 几个低次牛顿-科特斯求积公式的MATLAB实现先用M文件定义一个名为f1.m的函数： % i是要调用第几个被积函数g(i)，x是自变量 function f=f1(i,x) g(1)=sqrt(x); if x==0g(2)=1; elseg(2)=sin(x)/x; endg(3)=4/(1+x^2);f=g(i);程序一：function [C,g]=NCotes(a,b,n,m)% a，b分别为积分的上下限； % n是子区间的个数；% m是调用上面第几个被积函数；% 当n=1时计算梯形公式；当 i=n;h=(b-a)/i; z=0; for j=0:ix(j+1)=a+j*h; s=1; if j==0 s=s; else for k=1:j s=s*k; end end r=1;if i-j==0 r=r; elsefor k=1:(i-j) r=r*k; end endif mod((i-j),2)==1 q=-(i*s*r); elseq=i*s*r; end y=1;for k=0:i if k~=jy=y*(sym('t')-k); end endn=2时计算辛浦生公式，以此类推；l=int(y,0,i); C(j+1)=l/q;z=z+C(j+1)*f1(m,x(j+1)); endg=(b-a)*z1）当输入a?0,b?1，n?1,m?2时，即在MATLAB命令窗口输入 && NCotes(0,1,1,2)即可得用梯形公式的积分值和相应科特斯系数 如图3.12）当输入a?0,b?1，n?2,m?2时，即在MATLAB命令窗口输入 && NCotes(0,1,2,2)即可得用辛浦生公式的积分值和相应科特斯系数 如图3.23）当输入a?0,b?1，n?4,m?2时，即在MATLAB命令窗口输入 && NCotes(0,1,4,2)即可得用科特斯公式的积分值和相应科特斯系数 如图3.3图 3.1数值积分算法与MATLAB实现+毕业论文60_matlab数值积分图3.2图3.3第二节 复化求积公式的MATLAB实现一、复化梯形求积公式的MATLAB实现通过f(x)的n?1个等步长节点逼近积分?baf(x)dx?h2n?1(f(a)?f(b))?h?f(xk)k?1其中，xk?a?kh，x0?a,xn?b。程序二：function s=trapr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积；h=(b-a)/n;s=0;for k=1:(n-1)x=a+h*k;s=s+feval('f',x);endformat longs=h*(feval('f',a)+feval('f',b))/2+h*s;先用M文件定义一个名为f.m的函数：function y=f(x)if x==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入&& trapr1('f',0,1,4)回车得到 如图3.4图 3.4若取子区间的个数n?8在MATLAB命令窗口中输入&& trapr1('f',0,1,8)回车得到 如图3.5图3.5二、 复化辛浦生求积公式的MATLAB实现程序三：function s=simpr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积,即所求积分数值；h=(b-a)/(2*n);s1=0;s2=0;for k=1:nx=a+h*(2*k-1);s1=s1+feval('f',x);endfor k=1:(n-1)x=a+h*2*k;s2=s2+feval('f',x);ends=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3;先用M文件定义一个名为f.m的函数：function y=f(x)if x==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入&& simpr1('f',0,1,4)回车得到 如图3.6图3.6若取子区间个数n?8时在MATLAB命令窗口中输入&& simpr1('f',0,1,8)回车得到 如图3.7图3.7三、 复化科特斯求积公式的MATLAB实现程序四：function s=cotespr1(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积,即所求积分数值；h=(b-a)/n;s1=0;s2=0;s3=0;s4=0;for k=1:nx=a+(4*k-3)*h/4;s1=s1+feval('f',x);endfor k=1:n数值积分算法与MATLAB实现+毕业论文60_matlab数值积分x=a+(4*k-2)*h/4;s2=s2+feval('f',x);endfor k=1:nx=a+(4*k-1)*h/4;s3=s3+feval('f',x);endfor k=1:(n-1)x=a+4*k*h/4;s4=s4+feval('f',x);ends=h*(7*feval('f',a)+7*feval('f',b)+32*s1+12*s2+32*s3+14*s4)/90;同样先用M文件定义一个名为f.m的函数：function y=f(x)if x==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入&& cotespr1('f',0,1,4)回车得到 如图3.8图3.8第三节 龙贝格求积公式的MATLAB实现构造T数表来逼近积分?f(x)dx?R(J,J) ab其中。R(J,J)表示T数表的最后一行，最后一列的值。程序五：function [R,quad,err,h]=romber(f,a,b,n,delta)% f是被积函数% a,b分别是积分的上下限% n+1是T数表的列数% delta是允许误差% R是T数表% quad是所求积分值M=1;h=b-a;err=1J=0;R=zeros(4,4);R(1,1)=h*(feval('f',a)+feval('f',b))/2while ((err&delta)&(J&n))|(J&4)J=J+1;h=h/2;s=0;for p=1:Mx=a+h*(2*p-1);s=s+feval('f',x);endR(J+1,1)=R(J,1)/2+h*s;M=2*M;for K=1:JR(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);enderr=abs(R(J,J)-R(J+1,K+1));endquad=R(J+1,J+1)先用M文件定义一个名为f.m的函数：function y=f(x)if x==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入&& romber('f',0,1,5,0.5*(10^(-8)))回车得到 如图3.9图3.9第四节 高斯-勒让德求积公式的MATLAB实现程序六：function [A,x]=Guass1(N)i=N+1;f=((sym('t'))^2-1)^i;f=diff(f,i);t=solve(f);for j=1:ifor k=1:iX(j,k)=t(k)^(j-1);endif mod(j,2)==0B(j)=0;elseB(j)=2/j;endendX=inv(X);for j=1:iA(j)=0;x(j)=0;for k=1:iA(j)=A(j)+X(j,k)*B(k);x(j)=x(j)+t(j);endx(j)=x(j)/k;endfunction g= GuassLegendre (a,b,n,m)% a,b分别是积分的上下限；% n+1为节点个数；% m是调用f1.m中第几个被积函数；[A,x]=Guass1(n);g=0;for i=1:n+1y(i)=(b-a)/2*x(i)+(a+b)/2;f(i)=f1(m,y(i));g=g+(b-a)/2*f(i)*A(i);end用M文件分别把上面两个自定义函数定义为名为Guass1.m函数和GuassLegendre.m函数用M文件定义一个名为f1.m的函数function f=f1(i,x)g(1)=sqrt(x);if x==0g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x^2);f=g(i);在MATLAB命令窗口中输入&& GuassLegendre (0,1,2,2)&& GuassLegendre (0,1,3,2)回车得到 如图3.10图3.10第五节 各种求积算法的分析比较例1 分别用不同的方法计算积分I??10sinxx，并作比较。下面用各种求积公式分别计算积分，并给出了相应的计算误差，进行比较，结果如下：1、用Newton-Cotes公式当n=1时，即用梯形公式，用程序一数值积分算法与MATLAB实现+毕业论文60_matlab数值积分在MATLAB命令窗口中输入&& NCotes(0,1,1,2) 得I?0.95R?0.23当n=2时, 即用Simpson公式,用程序一在MATLAB命令窗口中输入&& NCotes(0,1,2,2)得I?0.59 R?0.407当n=4时, 即用科特斯公式，用程序一在MATLAB命令窗口中输入&& NCotes(0,1,4,2)得I?0.67 R?0.5132、用复化梯形公式令h=1/8=0.125，用程序二在MATLAB命令窗口中输入&& trapr1('f',0,1,8),得?10sinxx?h2?f(0)?2?f(h)???f(7h)??f(1)?=0.0R?0.1823、用复化辛浦生公式令h=1/8=0.125，用程序三在MATLAB命令窗口中输入&& simpr1('f',0,1,8),得?10sinxxdx?h3?f(0)?4?f(h)???f(7h)??2?f(2h)???f(6h)??f(1)?=0.95 R?0.7674、用Romberg公式 ，用程序五在MATLAB命令窗口中输入&& romber('f',0,1,5,0.5*(10^(-8)))，得?10sinxx?0.18 R?0.0025、用高斯-勒让德求积公式，令x?(t?1)/2，I??1?1sin(t?1)/2t?1（1） 用2个节点的Gauss公式I?0.33（2） 用3个节点的Gauss公式，用程序六在MATLAB命令窗口中输入&& GuassLegendre (0,1,2,2)，得I?0.473 R?0.290算法比较? 此例题的精确值为0.183? 由例题的各种求积算法可知：? 对Newton-cotes公式，当n=1时只有1位有效数字，当n=2时有3位有效数字，当n=4时有7位有效数字。? 对复化梯形公式有2位有效数字，对复化辛浦生公式有7位有效数字。 ? 用复化梯形公式，对积分区间[0,1]二分了11次用2049个函数值，才可得到7位准确数字。? 用Romberg公式对区间[0,1]二分3次，用了9个函数值，能得到同样的结果；二分4次，用了14个函数值，却能得到14位有效数字。? 用高斯-勒让德求积公式仅用了3个函数值，就能得到同样比较精确的6位有效数字。第六节 本章小结本章主要对各种求积算法，即牛顿-科特斯、复化求积、龙贝格、高斯-勒让德求积公式的算法在计算机上运用MATLAB软件进行编程实现，从而简化了运算。并且通过实例计算，误差分析对各种求积公式进行了分析比较。结 论本文主要讨论了数值积分的计算方法并通过MATLAB软件编程实现，通过前面的研究我们知道求数值积分近似值的计算方法很多，包括Newton-Cotes求积公式、复化求积公式、Romberg求积公式、高斯求积公式等等 。其中Newton-Cotes方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，Simpson方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；Cotes积分方法的误差比Simpson积分精度高两个数量级。因此，一般情况下，代数精度越高，积分公式计算精度也越高。但是高阶的Newton -Cotes方法的收敛性没有保证，因此，在实际计算中很少使用高阶的Newton -Cotes公式。复化梯形积分方法比单独的梯形积分精度高，它的积分精度和被积函数有关，还和复化积分时的步长有关。复化Simpson积分公式比单独的Simpson积分公式高近7个数量级，效果明显。Romberg方法收敛速度快、计算精度较高，但是计算量较大。Gauss求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦、而且要求已知积分函数。一般来说，Newton-Cotes方法的代数精度越高，数值积分的效果越好、越精确。当积分区间比较大的时候，可以采用复化积分方法可以得到更好的效果；变步长积分方法不仅可以很好地控制计算误差，并且可以寻找到适当的积分步长；Romberg积分方法可以更好地利用变步长复化积分公式得到的积分序列从而得到更为精确的数值结果，是一个较好的数值积分方法。高斯求积方法精确度高，收敛性快也是一种很优秀的数值积分方法。致 谢行文至此，我的这篇论文已接近尾声；岁月如梭，我四年的大学时光也即将敲响结束的钟声。离别在即，站在人生的又一个转折点上，心中难免思绪万千，一种感恩之情油然而生。首先感谢重庆邮电大学四年来对我的培养，是博学的老师们教会了我学习的方法、锻炼了我思考的能力、指明了我未来奋斗的方向，从而使我进一步明确了人生的目标。其次，我要感谢我的指导老师―郑继明教授，他的严谨细致、一丝不苟的作风一直是我工作、学习中的榜样；他的循循善诱的教导和不拘一格的思路给予我无尽的启迪。在撰写整个毕业论文的过程当中，他为我们考虑到了每一个细节，从开题报告到毕业论文的拟定修改上，郑老师更是不厌其烦的为我们做好每一步的细心指导。对此，我表示衷心地感谢。没有郑老师，我的论文也不可能这么顺利的完成。同时，我也要感谢每一位给过我帮助的老师和同学，在我撰写论文的过程当中同样给了我大量有益的建议，在此一并向他们表示真诚的感谢，感谢他们对我的支持和帮助。最后感谢这篇论文所涉及到的各位学者，本文引用了数位学者的研究文献，如果没有各位学者的研究成果带给我的的帮助和启发，我将很难完成本篇论文的写作。由于我的学术水平有限，所写论文难免有不足之处，恳请各位老师和学友批评指正。最后，衷心感谢评阅论文及参加答辩的各位老师！参考文献[1] 张德丰等．MATLAB数值计算方法［M］．北京：机械工业出版社，2010．[2] 郑继明等．数值计算方法［M］．重庆：重庆大学出版社，2005．[3] 樊守芳．Newton-Cotes数值求积公式的注记［J］． 枣庄学院学报，2011，28(2)：186-190．[4] 刘小伟．基于MATLAB的复合梯形数值积分法的研究与实验［J］． 甘肃联合大学学报(自然科学版)，)：20-23．[5] 余丹．用Simpson公式进行数值积分［J］． 华北电力大学数理系学报，2010，(23)：189-190．[6] Davis，P．J．and P．Rahinowitz．Methods of Numerical integration(secondedition)[M]． Academic Press. New York，1984．[7] 曹丽华． 一类广义Gauss型求积公式［J］． 深圳大学数学物理学报，2007，27(3)：524-534．[8] 郭晓斌．复化两点Gauss-Legendre公式及其误差分析［J］． 西北师范大学学报，)：49-51．[9] 王建强． 多种数值积分方法比较分析［J］． 武汉大学测绘轩辕学报，2010，2(1)：104-106．[10] 刘玉娟，陈应祖．龙贝格积分法及其应用编程［J］． 重庆科技学院学报，)：97-99．[11] 朱叶志．MATLAB数值分析与应用［M］．北京：机械工业出版社，2009．[12] Li F，Li X．The neighbor-scattering number can be computed in polynomialtime for interval graphs[J]．Computers and Mathematics with Applications，9-686．[13] 伍丽华，周玲丽．数学软件教程［M］．广州：中山大学出版社，2008．[14] 熊华，杨国孝．一类振荡函数的数值积分方法［J］． 北京理工大学学报，)：280-284．数值积分算法与MATLAB实现+毕业论文60_matlab数值积分[15] Varma．A K ．On optimal quadrature formulae［J］． Studia Sci Math Hungar，7-446．[16] 戴立辉．一些数值积分公式中间点的渐近性［J］． 闽江学院学报，2011，32(2)：5-7．[17] WANG X H，HE Y G，ZENG Z z．Numerical integration studybased ontriangle basis neural network algorithm[J]．Journal of Electronics&Information Technology．)：394-399．[18] GUO Ru．A Note on Newton-Cote?S Numerical Integral Formula［J］. Journalof Xinjiang University(Natural Science Edition)，)：186-190． [19] 刘小伟．基于MATLAB的变步长梯形数值积分法的研究与实验［J］．廊坊师范学院学报(自然科学版)，)：39-42.附 录一、英文原文：An Optimal Adaptative NumericalIntegration MethodA.Nicolet, A.Genon,W.Legros, M.Ume, F.DelThis author is a Research Assistant with the Belgian National Fund for Scientific Research University of Liege-Dept. of Electrical Engineering-Institut Montefiore Sart Tilman B28, B-4000 Liege,Belgium ABSTRACTThis paper presents an adaptative extension of the Gaussian integration method. It is well known that the Gaussian integration method is optimal for sufficiently smooth functions(i.e. which may be approximated by a polynomial) in the sense that it gives the maximum accuracy for a given number of nodes. Unfortunately it is not always possible to choose a priori the number of nodes for the integration. One alternative is to try successive Gaussian formulae with an increasing number of points until they agree with the required accuracy. In this case, most of the advantages of the method are lost. A less accurate but naturally adaptative method such as the Romberg method may become a better solution.The idea of the optimal adaptative method is to find a series of integration formulate with an increasing number of nodes in order that the set of abscissae of lower order formulae is a subset of abscissae of higher order formulae. Then, the sequential evaluation of formulae of increasing order only requires the addition of new points. Under this constraint, the remaining degrees of freedom(the newabscissa and all the weight factors)are used to obtain formulae of maximum order. INTRODUCTIONSome well-known and classical methods[1] are reviewed in order to situate the new method. Those methods are valid for sufficiently smooth integrands. Amongst the simplest, this method consists in choosing equally-spaced points between the endpoints and to approximate the function by piecewise linear functions. The trapezoidal rule is:I(h)??ba1?1?f(x)dx?h?f0?f1???fN?1?fN?2?2?(1)With: h?(b?a)/N and fi?f(a?ih) The basic idea is to use the results from successive refinements of the trapezoidal rule:I0?I(b?a),I1?I((b?a)/2),?,Ii?I((b?a)/2)? (2)iThe Richardson extrapolation is applied to this sequence in order to eliminate high order error terms[2].This method is naturally adaptative. For each refinement, a new trapezoidal approximation is computed with the number of points multiplied by two and reusing the previously computed values of the function. Then the Richardson extrapolation is applied to the new sequence. This process is repeated until the required accuracy is reached. The number of function evaluations is not known a priori and depends on the integrand. Most of the integration rules to integrate the following expression :?bak(x)f(x)dx (3)have the from :n?wi?1if(xi) (4)The approximation is a linear combination, with weight factors wi,of values of the function f(x), for n abscissae xi.In the rest of the paper, the integral (3), i.e. the integral on the interval [a,b] of the product of the function f(x) with the kernel k(x),will be referred as “the integral of f(x)”. Formulae (4) are tabulated in the literature for some kernelsk(x)and an associated interval [a,b].Orthogonal polynomials [3] associated to the set ?kernel k(x)-interval [a,b]? play an important role in the theory (table 1 ).Table 1. Orthogonal polynomials [4]The orthogonality of polynomialspi(x),is expressed by [3] :?baif i?j?o k(x)pi(x)pj(x)dx?? (5)? if i?j ?iA numerical integration rule is characterized by an integer p such that all thepolynomials Qk(x) of order k less than or equal to p are exactly integrated :?bank(x)Qk(x)dx??wQii?1k(xi) (6)Gaussian integration rules are optimal in the sense that p is maximal. An n point rule that has 2n degrees of freedom (n abscissae and n weight factors ) can integrate exactly all the polynomials up to the order 2n-1 (i.e. polynomials that have up to 2n coefficients).The necessary and sufficient condition is that all the powers of x up to the order 2n-1 are correctly integrated on the considered interval :?bank(x)xdx?mk?k?i?1wixi k?0,1?,kn,?2(7)The relations (7) constitute a system of 2n equations with 2n unknowns (thexiand the wi) whose solution gives the parameters of integration rule. Classicalapproach is to introduce an auxiliary polynomial ?(x) whose roots are the abscissae xi :nniii?(x)??(x?x)??cxi?1i?0(cn?1) (8)Sums of equations (7) weighted by the coefficients ciare constructed , in order to do appear ?(xi) :nn?i?0cimi?p??i?0?n?i?pci??wjxj???j?1?pj?ni?p?wj??cjxjj?1?i?0n?? (9)?n??w?(xj?1)xj?0 for p?0,?,n?1 jRelations (9) constitute a system of n equations with a unknowns that gives the coefficients ci :数值积分算法与MATLAB实现+毕业论文60_matlab数值积分n?1?cmii?0i?p?mn?p for p?0,?,n?1 (10 )The solution of the system (10) determines the polynomial?(x) whose roots are the abscissae xi of the integration rule.A relationship between this polynomial and the orthogonal polynomials associated to the problem may be found.pn(x)is the orthogonal polynomial of degree n associated to the problem.Any polynomial f2n?1(x) of degree 2n-1 may be expressed as :f2n?1(x)?g(x)pn(x)?r(x) (11 )With a quotient polynomialg(x) and a rest polynomial r(x), both of degree at most equal to n-1.The polynomial f2n?1(x) is integrated exactly :?bak(x)g(x)pn(x)dx??k(x)r(x)dx?abnniin?wg(x)pi?1(xi)??wir(xi) (12)i?1The polynomial pn(x) being orthogonal to all the polynomials of degree less or equal to n, the first term of the left hand member is equal to zero. The relation(12) is always true only if the first term of the right hand member is cancelled, what is assured if pn(xi) is equal to zero for i?1,?,n. Thus pn(x) is equal to ?(x) apart from a constant factor.The theory of the orthogonal polynomials assures that their roots are simple, real and situated in the interval [a,b]. The abscissae being determined, the wimust be computed in order that any polynomial of order less than or equal to n-1 is exactly integrated. This problem is solved in the next paragraph independently of the Gaussian method.If the number of points required to compute an integral with the satisfactoryaccuracy is known a priori, the Gaussian rule is the best method. Unfortunately ,this is not the case in most practical problems, A possibility is to try a sequence of Gaussian rules of increasing order until the difference between two approximations is less than the required accuracy. But, as the abscissae differ from one rule to the other, the advantages are lost and a naturally adaptative method such as the Romberg method becomes more relevant.WEIGHT COEFFICIENT DETERMINATIONThe problem of determining the ?optimal? weighting factors wi associated to any given n abscissae xi may be solved in a general way [5]. The n factors wi are computed in order that any polynomial up to order n-1 is integrated exactly. A polynomial fn?1(x) of order n-1 is integrated by the rule :?bank(x)fn?1(x)dx??wi?1ifn?1(xi) (13 )The Lagrange interpolation polynomials Li(x) are defined by :nLi(x)??j?1j?i(x?xj)(xi?xj)(14)Any polynomial fn?1(x) may be writtennfn?1(x)??L(x)fii?1n?1(xi) (15 )Introducing expression (15) in (13), gives, for any set of fn?1(xi) :?ba?n?k(x)??Li(x)fn?1(xi)?dx??i?1?n?wi?1ifn?1(xi) (16 )Identifying the coefficients of fn?1(xi) in the two members gives :wi??bak(x)Li(x)dx (17 )In the particular case of the Gaussian method, it can be shown that thisdetermination of the wi is equivalent to the solution of the n first equations of (7), linear with respect to wi, with the xi given. The wi obtained for the Gaussian rules are always positive if k(x) is positive on the interval [a,b]. Indeed, the polynomial L2(x) of degree 2(n-1) (the square of a Lagrange interpolating ipolynomial ) may be integrated exactly by the Gaussian rule :?bank(x)L(x)dx?2i?j?1Li(xj)wj (18)2The first member is positive become k(x)?0 on [a,b] and L2(x)?0. The iright hand member terms are all equal to zero except for i=j and than :wi??bak(x)Li(x)dxL(xi)2i2?0 (19 )This property assures a good response to Gaussian rules from the rounding error point of view.SINGULAR AND QUASI-SINGULAR INTEGRALSThe boundary element method involves the integration of singular and quasi-singular kernels. For the numerical integration method, the integrands were supposed to be smooth enough functions to be well approximated by polynomials.Nevertheless, the Patterson method may be efficiently applied if a change of variable is made in order to even out peaks or singularities[8,9].For instance, the integral of a function f(t) with a singularity or a peak att?tminis considered (in practice, this occurs when an influenced point is close to aninfluencing element, tmin corresponds to the parametric coordinate of the point of the influencing element at a minimum distance of or the same as the influenced point ) :I??10f(t)dt (20 )The following change of variable is performed :t?tmin?uThe integral (31) becomes :3(21)I?3?f(u?tmin)udu32(22 )In expression(22), the point corresponding to t?tminis u?0. Thus the singularity or the peak is eliminated by the term u2in (22).The practical algorithm is :- Check if the influenced point is on or close enough to the influencing element(The proximity criterion has been empirically chosen : a point is close enough to an element if its distance is less than one tenth of the length of the element ); - If the point is close enough, choose the expression(22), otherwise, chooseexpression(20);- Apply the adaptative Patterson method on the chosen expression. The criterionfor stopping the process is to have two successive approximations with a relative difference less than 10?4.CONCLUSIONOne of the characteristic of boundary element method is that it involves the computations of integrals ranging from very easy to singular. The proposed algorithm allows for the use of an efficient method with an adaptative order, and workable for the whole set of integrations. This provides accuracy and security(error is controlled) for a rather low computational cost.REFERENCES[1] P.J. DAVIS, P. RABINOWITZ. “Numerical integration”. Blaisdell,1967. [2] W.H. PRESS, B.P. FLANNERY, S.A. TEUKOLSKY, W.T. VETTERLING.“Numerical recipes : the art of scientific computing”. Cambridge University Press, 1989.[3] H. HOCHSTADT. “The functions of mathematical physics”. Dover, 1986. [4] M. ABRAMOWITZ, I.A.STEGUN ed.“Handbook of mathematical functions”.Dover, 1972.[5] T.N.L. PATTERSON. “The optimum addition of points to quadratureformulae”. Math. Comp. 22, 1968, pp.847-856.[6] T.N.L. PATTERSON. “On some Optimally and Lobatto based quadratureformulae”. Math. Comp. 22,1968, pp.877-881.[7] T.N.L. PATTERSON. “Algorithm for automatic numerical integration overa finite interval”. Communications of the ACM,Vol.16, no11,November 1973, pp.694-699.[8] J.P. ADRIAENS, P. BOURMANNE, F. DELINCE, A. GENON, A. NICOLET,W.LEGROS. “Numerical computations of eddy currents in thin plates”. IEEE Transactions on Magnetics,Vol.MAG-26,no5, September 1990, pp.. [9] A. NICOLET. “Modelisation du champ magnetique dans les systemscomprenant des milieu non lineaires”. PH.D. Thesis, University of Liege, 1991.数值积分算法与MATLAB实现+毕业论文60_matlab数值积分二、英文翻译：最优的适应性数值积分法A.Nicolet, A.Genon,W.Legros, M.Ume, F.Del作者是一个国家基金与比利时列日系的科研大学研究助理。电气工程学院 蒙蒂菲奥里SART蒂尔曼B28的B-4000列日，比利时 摘要本文提出了一种自适应高斯积分法的扩展。众所周知对于被积函数是充分光滑的函数（即可能由一个多项式近似代替），从某种意义上来说，高斯积分法是最佳的选择，它能给的最大精度为给定的大量节点数量个数。不幸的是，并不总是能选择为一体的节点数量。一种替代方法是尝试用不断增加的大量的节点，不断的连续用高斯公式，直到他们达到与所要求的精度。在这种情况下，该方法的大部分优点都将失去。一个不太准确，但自然适应性的方法，如龙贝格方法有可能成为一个更好的解决方案。最优自适应方法的思想是用大量不断增加的节点找到了一系列的整合公式，从而制定一套较低阶公式的横坐标是高阶公式的横坐标的一个子集。然后，为了增加公式的阶的序贯评价，只需要增加新的节点。在此约束下，剩余的自由度（新的横坐标和所有的权重因子）用于获得最大阶数的公式。引言对一些知名的经典方法[1]进行审查，以为了引出新方法。这些方法是有效的对于被积函数是充分光滑的函数。 梯形法梯形法是所有方法当中最简单的，这种方法选择同等的间隔点在端点和近似分段线性函数之间。梯形公式是 ：I(h)??ba1?1?f(x)dx?h?f0?f1???fN?1?fN? (1)2?2?其中 ：h?(b?a)/N 和 fi?f(a?ih) 龙贝格方法其基本思想是利用连续改进的梯形法的结果 ：I0?I(b?a),I1?I((b?a)/2),?,Ii?I((b?a)/2)? (2)i理查德森外推法应用到这个序列中，以消除高阶误差[2]。这种方法是自然适应的。对于每一个改进，每次节点数乘以二得到一个新的梯形法的近似计算结果，和重复使用先前计算的函数值。然后把理查德森外推法应用到新的序列中。反复进行这个过程，直至达到所需的精度。大量函数评估之前是不知道的这要取决于被积函数。一般规则大部分的整合规则如下列表达式：?bak(x)f(x)dx (3)有形式 ：n?wi?1if(xi) (4)这个近似的线性组合，权重因子wi， f(x)对于n个横坐标xi的函数值。 在本文的剩余部分，积分(3)，即被积函数f(x)和权函数k(x)在积分区间。公式(4)被制成表格对于一些[a,b]上的函数值的和，被称为“f(x)的积分”权函数k(x)和一个相关联的区间[a,b]，在这个文献中。正交多项式[3]相关集‘在区间[a,b]上的权函数k(x)’充当了一个重要的角色在这个理论中（表 1）。正交多项式[4] 表1.这个正交多项式pi(x)，被表述在[3]中?baif i?j?ok(x)pi(x)pj(x)dx?? (5)? if i?j ?i这种数值积分方法以所有小于或等于p的k阶多项式Qk(x)完全整合的整体p为特点：?bank(x)Qk(x)dx??wQii?1k(xi) (6)高斯方法从p是最大的意义上来说，高斯积分法是最优的选择。一个n点规则，有2n个自由度（n个横坐标和n个权重因子）可以完全地整合所有阶为2n-1的多项式（即多项式有2n个系数）。所有x的权限阶为2n-1次多项式能被正确地整合在这个被考虑的区间上的充分必要条件是：?bank(x)xdx?mk?k?i?1wixi k?0,1,?,2n?1 (7)k公式(7)关系构成一个具有2n个未知数2n个方程的体系(xi和wi)，其解决方法是提供集成方法的参数。经典的方法就是引入一个辅助多项式?(x)，其根是横坐标xi：nniii?(x)??(x?x)??cxi?1i?0(cn?1) (8)构造加权系数ci的方程(7)的总和，为了引出?(xi)：nn?i?0cimi?p??i?0?n?i?pci??wjxj???j?1?p?ni?p?wj??cjxjj?1?i?0n?? (9) ?n??j?1wj?(xj)xj?0 其中 p?0,?,n?1公式(9)构成一个具有系数ci的一个未知数的n个方程体系：n?1?cmii?0i?p?mn?p 其中 p?0,?,n?1 (10)这个体系(10)的解决方案是确定多项式?(x)的根是整合方法的横坐标xi 多项式与正交多项式之间的相关问题可以被发现。pn(x)是n阶相关问题的正交多项式。任何2n-1阶多项式f2n?1(x)可以被表述为：f2n?1(x)?g(x)pn(x)?r(x) (11)一个商数多项式g(x)和剩余多项式r(x)，都是最多能达到n-1阶。 多项式f2n?1(x)被完全地整合：?bak(x)g(x)pn(x)dx??k(x)r(x)dx?abnniin?wg(x)pi?1(xi)??wir(xi) (12)i?1多项式pn(x)与所有阶数小于或等于n的多项式正交，左手边成员的首项等于零。公式(12)总是正确的，即使右手边成员的首项被取消，如果对于i?1,?,n，pn(xi)被保证等于零，那么pn(x)等同于?(x)，等于一个常数因子。正交多项式的理论确保它们的根是简单的，真实的位于区间[a,b]上。横坐标被确定，这个wi因子必须被计算为了使任何阶小于或等于n-1的多项式完全被整合。高斯求积方法的独立性这个问题在下一段中将得到解决。如果知道区间中用于计算的大量的节点和需要一个完整的令人满意的精度，高斯求积方法是最好的方法。不幸的是，在大部分实际问题中并不是这样的，一种可能性是尝试用阶数不断增加的一系列高斯求积公式，直到两个近似值的差值是小于所要求的精度。但是，由于横坐标和其他有一个不同的规则，也失去了它的优势。一个自然适应性的方法，如龙贝格方法变得更加相关。 权系数的确定决定最优的任何n个横坐标xi的权因子wi的问题可以被解决用一种常规的方式[5]。这n个因子wi被计算为了使任何阶为n-1的多项式能被完全的整合。运用这种方法一个n-1次的多项式fn?1(x)被整合：?bank(x)fn?1(x)dx??wi?1ifn?1(xi) (13)拉格朗日插值多项式Li(x)被定义为：nLi(x)??j?1j?i(x?xj)(xi?xj) (14)任何多项式fn?1(x)可以被写成nfn?1(x)??L(x)fii?1n?1(xi) (15)把表达式(15)代入(13)中，使任何序列fn?1(xi)：?ba?n?k(x)??Li(x)fn?1(xi)?dx??i?1?n?wi?1ifn?1(xi) (16)数值积分算法与MATLAB实现+毕业论文60_matlab数值积分在给出的两组数中确定fn?1(xi)的系数wi??bak(x)Li(x)dx (17)在高斯求积方法的特殊情况下，与公式(7)中第一个方程的解相关的wi的确定值能被显示。考虑线性相关wi，用给定的xi。这个用高斯方法获得的wi总是确定的，如果k(x)在区间[a,b]上是确定的。实际上，2(n?1)阶多项式L2(x)i（一个拉格朗日插值多项式的平方）可以被完全的整合用高斯方法：?bank(x)L(x)dx?2i?Lj?12i(xj)wj (18)(x)?0。当这第一个数是确定的，因为在区间[a,b]上 k(x)?0 并且 L2ii?j时右边所有项总是等于零，然后：wi??bak(x)Li(x)dxLi(xi)22?0 (19)此性质从四舍五入误差的角度保证了高斯求积公式的良好效果。奇异和准奇异积分边界元法涉及的奇异和准奇异积分的整合。对于数值积分法，被积函数应该是足够光滑的函数能被多项式很好的代替。然而，帕特森方法可以被有效的应用，如果一个变量的变化是为了使顶点或奇点平坦 [8,9]。例如，函数f(t)用一个奇点或者在t?tmin的一个顶点求得的积分（在实践中，当一个受影响的点靠近一个影响元素出现， tmin相对应的参数坐标最小距离的影响元素的点或相同影响点）：I??10f(t)dt (20)接下来跟着变量的改变被执行：t?tmin?u (21) 3这个积分(31)变成：I?3?f(u?tmin)udu (22) 32在表达式(22)中，t?tmin相对应的点是u?0。因此这个奇点或者顶点被消除通过等式(22)中的u2。实际的算法是：- 检查影响点是否在或者足够接近影响元素（这个接近标准凭经验选择：一个点是足够接近元素，如果其距离小于这个元素的长度的十分之一）； - 如果这个点是足够接近的，选择表达式(22)，否则，选择表达式(20)； - 应用适应性帕特森法选择表达式。停止这一过程的标准是有两个连续的近似值相对差小于10?4。结论边界元法的特点之一是，它涉及的对于奇异积分的计算是很容易的。该算法允许使用适应性阶，以及一整套整合的可行有效的方法。这提供了准确性和安全性（误差控制）在一个相当低的计算成本下。参考文献[1] P.J. DAVIS, P. RABINOWITZ. “Numerical integration”. Blaisdell,1967.[2] W.H. PRESS, B.P. FLANNERY, S.A. TEUKOLSKY, W.T. VETTERLING.“Numerical recipes : the art of scientific computing”. Cambridge University Press, 1989.[3] H. HOCHSTADT. “The functions of mathematical physics”. Dover, 1986.[4] M. ABRAMOWITZ, I.A.STEGUN ed.“Handbook of mathematical functions”.Dover, 1972.[5] T.N.L. PATTERSON. “The optimum addition of points to quadratureformulae”. Math. Comp. 22, 1968, pp.847-856.[6] T.N.L. PATTERSON. “On some Optimally and Lobatto based quadratureformulae”. Math. Comp. 22,1968, pp.877-881.[7] T.N.L. PATTERSON. “Algorithm for automatic numerical integration overa finite interval”. Communications of the ACM,Vol.16, no11,November 1973, pp.694-699.[8] J.P. ADRIAENS, P. BOURMANNE, F. DELINCE, A. GENON, A. NICOLET,W.LEGROS. “Numerical computations of eddy currents in thin plates”. IEEE Transactions on Magnetics,Vol.MAG-26,no5, September 1990, pp..[9] A. NICOLET. “Modelisation du champ magnetique dans les systemscomprenant des milieu non lineaires”. PH.D. Thesis, University of Liege, 1991. 分享： >
“matlab数值积分”相关文章}