基础实验二 定积分数值计算
- 格式:doc
- 大小:103.43 KB
- 文档页数:3
第1篇一、实验目的1. 理解数值计算的基本概念和常用算法;2. 掌握Python编程语言进行数值计算的基本操作;3. 熟悉科学计算库NumPy和SciPy的使用;4. 分析算法的数值稳定性和误差分析。
二、实验内容1. 实验环境操作系统:Windows 10编程语言:Python 3.8科学计算库:NumPy 1.19.2,SciPy 1.5.02. 实验步骤(1)Python编程基础1)变量与数据类型2)运算符与表达式3)控制流4)函数与模块(2)NumPy库1)数组的创建与操作2)数组运算3)矩阵运算(3)SciPy库1)求解线性方程组2)插值与拟合3)数值积分(4)误差分析1)舍入误差2)截断误差3)数值稳定性三、实验结果与分析1. 实验一:Python编程基础(1)变量与数据类型通过实验,掌握了Python中变量与数据类型的定义方法,包括整数、浮点数、字符串、列表、元组、字典和集合等。
(2)运算符与表达式实验验证了Python中的算术运算、关系运算、逻辑运算等运算符,并学习了如何使用表达式进行计算。
(3)控制流实验学习了if-else、for、while等控制流语句,掌握了条件判断、循环控制等编程技巧。
(4)函数与模块实验介绍了Python中函数的定义、调用、参数传递和返回值,并学习了如何使用模块进行代码复用。
2. 实验二:NumPy库(1)数组的创建与操作通过实验,掌握了NumPy数组的基本操作,包括创建数组、索引、切片、排序等。
(2)数组运算实验验证了NumPy数组在数学运算方面的优势,包括加、减、乘、除、幂运算等。
(3)矩阵运算实验学习了NumPy中矩阵的创建、操作和运算,包括矩阵乘法、求逆、行列式等。
3. 实验三:SciPy库(1)求解线性方程组实验使用了SciPy库中的线性代数模块,通过高斯消元法、LU分解等方法求解线性方程组。
(2)插值与拟合实验使用了SciPy库中的插值和拟合模块,实现了对数据的插值和拟合,并分析了拟合效果。
第一讲定积分的数值计算【主要目的】围绕定积分的概念与数值计算方法这一大家非常熟悉的主题,突出数值实验、几何观察、数值分析等实验特性,学生通过实验与理论的对照,加深对数学思想和数学知识的理解和掌握,学习如何从实验角度创新知识、发现知识,并上升到理论分析的角度。
【主要内容】定积分的数值计算方法,包括:矩形法、梯形法与辛普森法;对误差的了解:精度与收敛速度引言首先回忆一下函数在区间上的定积分概念的建立过程。
考虑在区间内任意插入个分点的分法:把分割成个小区间,第个子区间的长度为;任取数,做乘积,把所有这些乘积相加得到和式.如果无论区间怎样划分及分点怎样选取,当时,该和式都趋于同一常数,则称函数在区间上可积,且称此常数为在区间上的定积分,即。
称和式为积分和或黎曼和。
在定积分的概念中包含了两个任意性,即对区间的分割和点的选取都是任意的。
显然,对于区间的不同分割或者点的选取不同,得到的和式一般不同。
定积分的定义中要求在对区间无限细分()的条件下,所有这些和式都趋于同一数值。
这一点初学者较难理解。
我们将通过数值实验来加以理解。
当在区间上连续,为在区间上的原函数时,我们可以用牛顿-莱布尼兹公式方便地求得。
但是有些函数其原函数不能用初等函数表示出来,这样对应的定积分通常也不能用牛顿-莱布尼兹公式算出其精确值。
而且,在自然科学与工程技术中有许多问题,被积函数并不是用具体函数表达式解析表示的,而经常是通过实验或测量方法用表格或图形给出的,这就导出了定积分的数值计算问题。
我们将利用“分割取近似,作和求极限”这一定积分思想方法,来构造一些数值计算方法,并进行数值实验。
实验一定积分概念的深化——达布和设函数在区间上有界。
考虑将将区间任意分割成个子区间()的分法,设在子区间上的上、下确界分别为,称为在子区间上的振幅,和式分别称为关于该分割的达布(Darboux)大和与达布小和。
由定义可知,函数对应于同一分割的积分和有无穷多个,但达布大和与达布小和却都各只有一个。
几种定积分的数值计算方法一、梯形法则(Trapezoidal Rule):梯形法则是一种常见的确定积分的数值计算方法。
它的基本思想是通过将函数曲线上的曲线段看作是一系列梯形,然后计算这些梯形的面积之和来近似表示定积分的值。
具体来说,我们将定积分区间[a,b]均匀地划分为n个小区间,每个小区间的宽度为h=(b-a)/n,然后计算每个小区间内的梯形面积,再将这些面积相加即可得到定积分的近似值。
梯形法则的公式如下:∫(a to b) f(x) dx ≈ h/2 * (f(a) + 2f(a+h) + 2f(a+2h) + ... + 2f(a+(n-1)h) + f(b))梯形法则的优点是简单易懂,容易实现,并且对于一般的函数都能达到较好的近似效果。
然而,它的缺点是精度较低,需要较大的划分数n才能得到较准确的结果。
二、辛普森法则(Simpson's Rule):辛普森法则是一种比梯形法则更高级的确定积分方法,它通过将函数曲线上的曲线段看作是由一系列抛物线组成的,然后计算这些抛物线的面积之和来近似表示定积分的值。
与梯形法则类似,我们将定积分区间[a,b]均匀地划分为n个小区间,每个小区间的宽度为h=(b-a)/n,然后计算每两个相邻小区间内的抛物线面积,再将这些面积相加即可得到定积分的近似值。
辛普森法则的公式如下:∫(a to b) f(x) dx ≈ h/3 * (f(a) + 4f(a+h) + 2f(a+2h) +4f(a+3h) + ... + 2f(a+(n-2)h) + 4f(a+(n-1)h) + f(b))辛普森法则相较于梯形法则具有更高的精度,尤其对于二次或更低次的多项式函数来说,可以得到非常准确的结果。
但是,辛普森法则在处理高次多项式或非多项式函数时可能会出现误差较大的情况。
三、高斯求积法(Gaussian Quadrature):高斯求积法是一种基于插值多项式的数值积分方法。
定积分的计算定积分是微积分中的一个重要概念,用来计算曲线与x轴之间的面积或曲线的弧长等问题。
本文将介绍定积分的概念、性质和计算方法。
一、定积分的概念定积分是一种数学运算,用来计算曲线与x轴之间的面积。
它的定义是在一个区间上划分出无穷多个小矩形,然后计算这些小矩形的面积之和,然后取极限。
用符号表示为∫f(x)dx,其中f(x)是被积函数,dx表示微元长度。
二、定积分的性质1. 定积分具有可加性,即∫(f(x) + g(x))dx = ∫f(x)dx + ∫g(x)dx。
2. 定积分的区间可加性,即∫[a,b]f(x)dx = ∫[a,c]f(x)dx + ∫[c,b]f(x)dx。
3. 定积分的值与被积函数的符号无关,即∫[a,b]f(x)dx = -∫[b,a]f(x)dx。
4. 定积分的值与积分区间的长度无关,即∫[a,b]f(x)dx = ∫[ka,kb]f(x)dx,其中k为任意非零常数。
三、定积分的计算方法计算定积分的方法有很多种,以下是一些常用的方法:1. 几何方法:对于一些简单的几何图形,我们可以利用几何的知识来求解。
例如,对于一个矩形的面积,可以直接计算长度乘以宽度。
2. 切割方法:将区间切割成无穷多个小区间,并计算每个小区间上的面积之和。
当小区间趋近于无穷小时,这个和就是定积分的值。
这种方法也被称为黎曼和的定义。
3. 牛顿-莱布尼茨公式:若函数G(x)是f(x)的一个原函数,则定积分可以通过G(b) - G(a)来计算,其中a、b是积分区间的端点。
4. 变量代换法:对于一些复杂的函数,可以通过变量代换来简化问题。
例如,对于∫(x^2 + 1)dx,我们可以令u = x^2 + 1,然后计算∫udu,最后再带回原来的变量。
5. 分部积分法:对于一些产品的积分,可以利用分部积分公式来求解。
该公式为∫u(x)v'(x)dx = u(x)v(x) - ∫v(x)u'(x)dx。
一、实验目的1. 理解积分的概念和基本性质。
2. 掌握数值积分的方法,包括矩形法、梯形法、辛普森法等。
3. 通过实际计算,加深对积分概念的理解。
二、实验原理积分是微积分学中的一个基本概念,表示一个函数在某区间内的累积变化量。
数值积分是指利用数值方法求解积分,常见的方法有矩形法、梯形法、辛普森法等。
1. 矩形法:将积分区间分成若干等份,用每个小区间的宽度乘以函数在该区间的值,再将所有小区间的乘积相加,得到积分的近似值。
2. 梯形法:将积分区间分成若干等份,用每个小区间的宽度乘以函数在该区间的平均值,再将所有小区间的乘积相加,得到积分的近似值。
3. 辛普森法:将积分区间分成若干等份,用每个小区间的宽度乘以函数在该区间的二次多项式近似值,再将所有小区间的乘积相加,得到积分的近似值。
三、实验步骤1. 选择一个具体的积分问题,例如:计算函数f(x) = x^2在区间[0,1]上的积分。
2. 根据所选择的积分方法,设置相应的参数。
例如,对于矩形法,需要设置小区间的数量n;对于梯形法,需要设置小区间的数量n;对于辛普森法,需要设置小区间的数量n。
3. 计算每个小区间的宽度,例如,对于区间[0,1],小区间的宽度为h = (1-0)/n。
4. 根据所选的积分方法,计算积分的近似值。
5. 比较不同积分方法的近似值,分析误差来源。
四、实验结果与分析以函数f(x) = x^2在区间[0,1]上的积分为例,进行数值积分实验。
1. 矩形法:取n=4,计算得到积分的近似值为0.5625。
2. 梯形法:取n=4,计算得到积分的近似值为0.6667。
3. 辛普森法:取n=4,计算得到积分的近似值为0.6667。
通过比较不同积分方法的近似值,可以发现辛普森法的误差较小,且随着n的增大,误差逐渐减小。
这表明辛普森法在数值积分中具有较高的精度。
五、实验总结1. 本实验通过数值积分方法,计算了函数f(x) = x^2在区间[0,1]上的积分,加深了对积分概念的理解。
几种定积分的数值计算方法数值计算定积分是计算定积分的一种近似方法,适用于无法通过代数方法求得精确解的定积分。
本文将介绍几种常见的数值计算定积分的方法。
1.矩形法(矩形逼近法):矩形法是最简单的数值计算定积分方法之一、它将定积分区间划分为若干个小区间,然后在每个小区间上取一个样本点,将每个小区间上的函数值乘以小区间的宽度,得到小矩形的面积,最后将这些小矩形的面积相加即可得到定积分的近似值。
矩形法有两种主要的实现方式:左矩形法和右矩形法。
左矩形法使用每个小区间的左端点作为样本点,右矩形法则使用右端点。
2.梯形法(梯形逼近法):梯形法是另一种常见的数值计算定积分方法。
它将定积分区间划分为若干个小区间,然后在每个小区间上取两个样本点,分别作为小区间的端点。
接下来,计算每个小区间上的函数值,然后将每个小区间上的函数值与两个端点连线所构成的梯形的面积相加,得到所有梯形的面积之和,最后得到近似的定积分值。
3.辛普森法:辛普森法是一种更为精确的数值计算定积分方法。
它将定积分区间分为若干个小区间,然后用二次多项式逼近每个小区间上的函数曲线。
在每个小区间上,辛普森法使用三个样本点,将函数曲线近似为一个二次多项式。
然后,对于每个小区间,计算该二次多项式所对应的曲线下梯形区域的面积,并将所有小区间的面积相加,得到近似的定积分值。
4. 龙贝格法(Romberg integration):龙贝格法是一种迭代的数值计算定积分方法,通过进行多次计算,逐步提高近似的精确度。
龙贝格法首先使用梯形法或者辛普森法计算一个初始近似值,然后通过迭代的方式进行优化。
在每次迭代中,龙贝格法先将区间划分成更多的子区间,并在每个子区间上进行梯形法或者辛普森法的计算。
然后,利用这些计算结果进行Richardson外推,从而得到更精确的定积分近似值。
通过多次迭代,龙贝格法可以逐步提高逼近的精确度。
上述介绍的四种数值计算定积分的方法都有各自的优势和适用范围。
实验二定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值,但它仅适用于被积函数的原函数能用初等函数表达出来的情形.如果这点办不到或者不容易办到,这就有必要考虑近似计算的方法.在定积分的很多应用问题中,被积函数甚至没有解析表达式,可能只是一条实验记录曲线,或者是一组离散的采样值,这时只能应用近似方法去计算相应的定积分.本实验将主要研究定积分的三种近似计算算法:矩形法、梯形法、抛物线法.对于定积分的近似数值计算,Matlab有专门函数可用.二、相关函数(命令)及简介1.sum(a):求数组a的和.2.format long:长格式,即屏幕显示15位有效数字.(注:由于本实验要比较近似解法和精确求解间的误差,需要更高的精度).3.double():若输入的是字符则转化为相应的ASCII码;若输入的是整型数值则转化为相应的实型数值.4.quad():抛物线法求数值积分.格式: quad(fun,a,b) ,注意此处的fun是函数,并且为数值形式的,所以使用*、/、^等运算时要在其前加上小数点,即 .*、./、.^等.例:Q = quad('1./(x.^3-2*x-5)',0,2);5.trapz():梯形法求数值积分.格式:trapz(x,y)其中x为带有步长的积分区间;y为数值形式的运算(相当于上面介绍的函数fun)例:计算x=0:pi/100:pi;y=sin(x);trapz(x,y)6.dblquad():抛物线法求二重数值积分.格式:dblquad(fun,xmin,xmax,ymin,ymax),fun可以用inline定义,也可以通过某个函数文件的句柄传递.例1:Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2,通过计算,比较Q1 与Q2结果(或加上手工验算),找出积分变量x、y的上下限的函数代入方法.Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi) 例2:Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件integrnd.m:function z = integrnd(x, y)z = y*sin(x);7.fprintf(文件地址,格式,写入的变量):把数据写入指定文件.例:x = 0:.1:1;y = [x; exp(x)];fid = fopen('exp.txt','w'); %打开文件fprintf(fid,'%6.2f %12.8f\n',y); %写入fclose(fid) %关闭文件8.syms 变量1 变量2 …:定义变量为符号.9.sym('表达式'):将表达式定义为符号.解释:Matlab中的符号运算事实上是借用了Maple的软件包,所以当在Matlab中要对符号进行运算时,必须先把要用到的变量定义为符号.10.int(f,v,a,b):求f关于v积分,积分区间由a到b.11.subs(f,'x',a):将 a 的值赋给符号表达式 f 中的 x,并计算出值.若简单地使用subs(f),则将f的所有符号变量用可能的数值代入,并计算出值。
第1篇一、实验目的1. 理解定积分的概念,掌握定积分的计算方法。
2. 熟悉数值积分的方法,提高数值计算能力。
3. 通过实验,验证定积分的计算结果,加深对定积分理论的理解。
二、实验原理定积分是数学分析中的一个基本概念,它表示函数在某一区间上的累积效果。
对于给定的函数f(x),在区间[a, b]上的定积分可以表示为:∫[a, b] f(x) dx其中,dx表示无穷小的区间宽度。
在实际计算中,定积分往往采用数值积分的方法进行近似计算。
三、实验仪器与软件1. 仪器:计算机2. 软件:MATLAB四、实验步骤1. 输入函数表达式:在MATLAB中输入待积分函数的表达式,例如f(x) = x^2。
2. 设置积分区间:设定积分的上下限a和b。
3. 选择数值积分方法:MATLAB提供了多种数值积分方法,如梯形法、辛普森法、高斯法等。
根据需要选择合适的方法。
4. 进行数值积分计算:调用MATLAB的数值积分函数,如quad函数,进行积分计算。
5. 结果分析:观察计算结果,与理论值进行对比,分析误差来源。
五、实验数据及结果1. 函数表达式:f(x) = x^22. 积分区间:[0, 1]3. 数值积分方法:辛普森法4. 计算结果:I ≈ 1.1666666667六、误差分析1. 理论值:∫[0, 1] x^2 dx = [x^3/3] |[0, 1] = 1/32. 误差来源:a. 数值积分方法的误差:由于数值积分方法是一种近似计算方法,其计算结果与真实值存在一定的误差。
b. 计算过程中的舍入误差:在计算过程中,由于计算机的浮点数表示,可能导致舍入误差。
3. 误差分析:计算结果与理论值相差较大,说明数值积分方法的误差较大。
在实际应用中,可以根据需要选择合适的数值积分方法,以减小误差。
七、实验结论1. 通过本次实验,掌握了定积分的计算方法,了解了数值积分的方法及其优缺点。
2. 了解了数值积分方法在计算过程中的误差来源,为实际应用提供了参考。
一、实验目的1. 熟悉数值计算的基本原理和方法。
2. 掌握常用的数值计算算法及其应用。
3. 提高数值计算软件的使用能力。
4. 培养分析问题和解决问题的能力。
二、实验环境1. 操作系统:Windows 102. 编程语言:Python3. 数值计算软件:NumPy、SciPy、Matplotlib三、实验内容1. 实验一:数值积分(1)实验目的:学习数值积分方法,计算定积分的近似值。
(2)实验内容:a. 使用辛普森法则计算函数f(x) = x^2在区间[0, 1]上的定积分。
b. 使用梯形法则计算函数f(x) = e^x在区间[0, 1]上的定积分。
(3)实验步骤:a. 编写Python代码,实现辛普森法则和梯形法则。
b. 分别使用两种方法计算定积分的近似值。
c. 对比两种方法的计算结果,分析误差来源。
2. 实验二:数值微分(1)实验目的:学习数值微分方法,计算函数在某点的导数近似值。
(2)实验内容:a. 使用中心差分法计算函数f(x) = sin(x)在x = π/2处的导数近似值。
b. 使用前向差分法和后向差分法计算函数f(x) = cos(x)在x = 0处的导数近似值。
(3)实验步骤:a. 编写Python代码,实现中心差分法、前向差分法和后向差分法。
b. 分别使用三种方法计算导数的近似值。
c. 对比三种方法的计算结果,分析误差来源。
3. 实验三:线性方程组求解(1)实验目的:学习线性方程组求解方法,掌握高斯消元法和迭代法。
(2)实验内容:a. 使用高斯消元法求解线性方程组:3x + 2y - z = 72x - y + 3z = -1-x + 2y + 2z = 4b. 使用雅可比迭代法求解线性方程组:3x + 2y - z = 72x - y + 3z = -1-x + 2y + 2z = 4(3)实验步骤:a. 编写Python代码,实现高斯消元法和雅可比迭代法。
b. 分别使用两种方法求解线性方程组。
第1篇一、实验目的1. 理解定积分的概念和意义。
2. 掌握定积分的计算方法。
3. 通过实验,加深对定积分理论知识的理解和应用。
二、实验原理定积分是微积分学中的一个基本概念,它表示某一曲线与x轴及两条垂直于x轴的直线所围成的平面图形的面积。
定积分的计算方法主要有直接积分法、分部积分法、换元积分法等。
三、实验仪器与材料1. 计算机2. 计算软件(如MATLAB、Mathematica等)3. 数学教材和参考书籍四、实验内容1. 理论学习(1)复习定积分的定义、性质和计算方法。
(2)了解定积分在物理学、工程学、经济学等领域的应用。
2. 实验操作(1)计算以下定积分:∫(0到1)x^2 dx(2)使用MATLAB软件进行定积分计算,并比较与手工计算的结果。
3. 结果分析(1)分析手工计算与MATLAB计算的结果差异,找出原因。
(2)总结定积分计算方法在实际应用中的优缺点。
五、实验步骤1. 理论学习(1)阅读数学教材,了解定积分的定义、性质和计算方法。
(2)查阅相关资料,了解定积分在各个领域的应用。
2. 实验操作(1)使用计算器或手工计算定积分∫(0到1)x^2 dx,得到结果。
(2)打开MATLAB软件,编写以下代码进行定积分计算:```matlabsyms x;I = int(x^2, 0, 1);```运行代码,得到结果。
3. 结果分析(1)比较手工计算和MATLAB计算的结果,发现两者基本一致。
(2)分析定积分计算方法在实际应用中的优缺点,总结如下:优点:- 定积分计算方法可以处理各种复杂函数,具有广泛的适用性。
- 定积分在物理学、工程学、经济学等领域的应用非常广泛。
缺点:- 对于某些复杂函数,计算过程较为繁琐,需要花费较多时间。
- 在实际应用中,可能需要根据具体情况选择合适的计算方法。
六、实验总结通过本次实验,我们了解了定积分的概念、性质和计算方法,并掌握了使用MATLAB软件进行定积分计算的方法。
基础实验二 定积分数值计算一、实验目的学习定积分的数值计算方法,理解定积分的定义,掌握牛顿-莱布尼兹公式。
二、实验材料2.1定积分的数值计算计算定积分⎰b adx x f )(的近似值,可将积分区间n 等分而得矩形公式 n a b n a b i a f dx x f n i b a ---+≈∑⎰=])1([)(1或n ab n a b i a f dx x f n i b a --+≈∑⎰=][)(1 也可用梯形公式近似计算na b b f a f n a b i a f dx x f n i ba -++-+≈∑⎰-=]2)()()([)(11 如果要准确些,可用辛普森公式n ab b f a f a b i a f n a b i a f dx x f n i n i b a 6)]()()2)21((4)(2[)(111-++--++-+≈∑∑⎰=-= 对于⎰10sin xdx ,矩形公式、梯形公式、辛普森公式的Mathematica 程序为a=0;b=1;k=10; f[x_]:=Sin[x];d=N[Integrate[f[x],{x,a,b}],k];(计算精确值)s1[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,0,m-1}],k];(取小区间左端点的矩形公式) s2[m_]:=N[Sum[f[a+(i+1/2)*(b-a)/m]*(b-a)/m,{i,0,m-1}],k]; (取小区间中点的矩形公式)s3[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,1,m}],k]; (取小区间右端点的矩形公式) s4[m_]:=N[Sum[(f[a+i*(b-a)/m]+f[a+(i+1)*(b-a)/m])/2*(b-a)/m,{i,0,m-1}],k]; (梯形公式)s5[m_]:=N[(b-a)/m/6*((f[a]+f[b])+2*Sum[f[a+i*(b-a)/m],{i,1,m-1}]+4*Sum[f[a+(i-1/2)*(b-a)/m],{i,1,m}]),k];(辛普森公式)r1[m_]:=d-s1[m];r2[m_]:=d-s2[m];r3[m_]:=d-s3[m];r4[m_]:=d-s4[m];r5[m_]:=d-s5[m];(误差)t=Table[{s1[m],r1[m],s2[m],r2[m],s3[m],r3[m],s4[m],r4[m],s5[m],r5[m]},{m,100,1000,100}]利用以上程序计算⎰10dx 、⎰10xdx 、⎰102dx x ,并对几个公式比较。
毕业论文文献综述信息与计算科学定积分的数值计算方法一、 前言部分在科学与工程计算中,经常要计算定积分()()().baI f f x dx a b =-∞≤≤≤∞⎰ (1.1)这个积分的计算似乎很简单,只要求出f 的原函数F 就可以得出积分(1.1)的值,即()()().I f F b F a =- (1.2)如果原函数F 非常简单又便于使用,那么式(1.2)就提供了计算起来最快的积分法.但是,积分过程往往将导出新的超越函数,例如,简单积分1dx x ⎰可引出对数函数,它已不是代数函数了;而积分2x edx -⎰,将引出一个无法用有限个代数运算、对数运算或指数运算组合表示的函数.有些积分虽然容易求解,并且原函数仍然是一个初等函数,但可能过于复杂,以致于人们采用(1.2)来计算之前还得三思而行[1].例如411dx C x =++⎰, (1.3) 采用式(1.3)这种“精确”表达式时,所需运算次数是个根本问题.由式(1.3)看出,需计算对数和反正切,因此只能计算到一定的近似程度.因此可以看出,这类表面上是“精确”的方法,实际上也是近似的.因此,我们常常需要探讨一些近似计算定积分的数值方法[2].通过人们的研究和发现,得出了很多数值计算的方法,比如利用牛顿-科茨求积公式,复合求积公式,龙贝格积分法,高斯求积公式,切比雪夫求积法等来解决定积分的数值计算问题.构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数,由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-柯茨公式,例如梯形公式与抛物线公式就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中,对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法,它具有公式简练、计算结果准确、使用方便、稳定性好等优点,因此在等距情形宜采用龙贝格求积公式.当用不等距节点进行计算时,常用高斯型求积公式计算,它在节点数目相同情况下,准确程度较高,稳定性好,而且还可以计算无穷积分[3].二、 主题部分2.1 牛顿-科茨求积公式[4]2.1.1 公式的一般形式[4]将积分(1.1)中的积分区间[],a b 分成n 等分,其节点k x 为1,()k x a kh h b a n=+=- (0,1,,)k n =L . 对于给定的函数f ,在节点k x (0,1,,)k n =L 上的值()k f x 为已知.那么f 在n+1个节点01,,,n x x x L 上的n 次代数插值多项式为00()().n nj n kk j k j j k x x p x f x x x ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ 如果记x a th =+,则上式可以写为00()().n nn kk j j k t j p x f x k j ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ (2.1) 在积分(1.1)中的被积函数f 用其n+1个节点的代数插值多项式()n p x 来代替,可 得 ()()()()bbn n aaI f f x dx I f p x dx =≈=⎰⎰.多项式的积分是容易求出的,因此把上式写为()()()nn n k k I f I f A f x =≈=∑, (2.2)其中 ()00(),n n n k k j j kb a t j A dt b ac n k j=≠--==--∏⎰ (2.3) ()00(1)().!()!n kn n n kj j kct j dt k n k n -=≠-=--∏⎰ (2.4) 公式(2.2)称为牛顿-科茨求积公式或称为等距节点求积公式,k A 称为求积公式系数,()n k c 称为科茨求积系数.牛顿-科茨求积公式的误差估计()n E f ()()n I f I f =-,由下面定理给出 定理2.1 (1) 如果n 为偶数,(2)n f +在[],a b 上连续,则有[]3(2)()(),,n n n n E f c hf a b ηη++=∈, (2.5)其中 201(1)(2)()(2)!n n c t t t t n dt n =---+⎰L . (2) 如果n 为奇数,(1)n f+在[],a b 上连续,则有[]2(1)()(),,n n n n E f c h f a b ηη++=∈, (2.6)其中 01(1)(2)()(1)!n n c t t t t n dt n =---+⎰L . 定义2.1 如果求积公式()()nbk k ak f x dx A f x =≈∑⎰对所有次数不高于n 的代数多项式等式精确成立,但存在n+1次的代数多项式使等式不成立,则称上式求积公式具有n 次代数精度.由定理2.1可知,牛顿-科茨求积公式(2.2)的代数精度至少是n 次,而当n 是偶数时,(2.2)的代数精度可达n+1次.2.1.2 梯形公式[5]在牛顿-科茨公式(2.2)中,取n=1时(1)(1)011,2c c ==所以有 []1()()()().2b aI f I f f a f b -≈=+ (2.7) 公式(2.7)称为梯形公式,如果用连接(),()a f a 和(),()b f b 的直线来逼近f ,并对这线性函数进行积分可得到1()I f .再用1()I f 来逼近()I f . 定理 2.2 若[]2,f Ca b ∈,则梯形公式(2.7)的误差为[]3111()()()()''(),,.12E f I f I f b a f a b ηη=-=--∈ 2.1.3 辛普森公式[6]在牛顿-科茨公式(2.2)中,取n=2,则有220011(1)(2),46c t t dt =--=⎰221014(2),26c t t dt =--=⎰ 222011(1),46c t t dt =-=⎰有此得到2()()()4()().32h a b I f I f f a f f b +⎡⎤≈=++⎢⎥⎣⎦(2.8) 其中1()2h b a =-.式(2.8)称为辛普森公式. 定理2.3 若[]4,f Ca b ∈,则辛普森公式(2.8)的误差为[]5(4)221()()()(),,.90E f I f I f h f a b ηη=-=-∈2.2 复化求积公式[7]上面已经给出了计算积分()()baI f f x dx =⎰的3个基本的求积公式:梯形公式,辛普森公式,牛顿-科茨公式,并给出了它们误差的表达式.由这些表达式可知其截断误差依赖于求积区间的长度.若积分区间的长度是小量的话,则这些求积公式的截断误差是该长度的高阶小量.但若积分区间的长度比较大,直接使用这些公式,则精度难以保证.为了提高计算积分的精度,可把积分区间分为若干个小区间,()I f 等于这些小区间上的积分和,然后对每个小区间上的积分应用上述求积公式,并把每个小区间上的结果累加,所得到的求积公式称为复化求积公式.将积分区间[],a b 作n 等分,并记,,0,1,,k b ah x a kh k n n-==+=L ,于是 11()()k kn x x k I f f x dx +-==∑⎰.2.2.1 复化梯形求积公式[8]如果需要求出一个已知函数()f x 在一个很大区间[],a b 上的积分,那么我们可以把区间分成n 个长度为x h ∆=的小区间,对每一个小区间用梯形法则,然后再把这些小区间上的积分值相加.于是就得到了计算定积分的复化梯形公式:1101210()()(222)22n bi i n n ai h hf x dx f f f f f f f -+-=≈+=+++++∑⎰L (2.9)整体积分误差等于n 个小区间上的积分误差之和:整体误差= []312''()''()''()12n h f f f ξξξ-+++L ,其中i ξ是第i 个小区间上的某一点.如果''()f x 在区间[],a b 上连续,那么由连续函数的性质可知,在区间[],a b 上存在点ξ使得''()i f ξ的平均值等于()f ξ.于是由于nh b a =-,有整体误差= 322''()''()()1212nh b a f h f O h ξξ--=-=, 局部误差是3()O h ,整体误差是2()O h .2.2.2 复化辛普森求积公式[9]对于积分()baf x dx ⎰,将[],a b 等分,每个小区间长度b ah n-=,节点记为 (0,1,2,,)k x a kh k n =+=L ,第k 个小区间记为[]1,(1,2,,)k k x x k n -=L .记[]1,k k x x -的中点为1121()2k k k xx x --=+,则复化辛普森公式为 1112()()()4()()6n bk k ak k h f x dx S h f x f x f x --=⎡⎤≈=++⎢⎥⎣⎦∑⎰.2.3 龙贝格积分[10]现在要介绍用龙贝格(Romberg )命名的一个算法,龙贝格首先给出了这种算法的递推形式,假设需要积分()baI f x dx =⎰ (2.10)的近似值.在讨论过程中函数()f x 和区间[],a b 将保持不变.2.3.1 递推梯形法则[10]设()T n 表示在长度是()/h b a n =-的n 个子区间上积分I 的梯形法则.根据()''()nbai f x dx h f a ih =≈+∑⎰,我们有 00()()''()''()nn n i i b a b a T h f a ih f a i n n ==--=+=+∑∑, (2.11) 这里求和符号中的两撇表示和式中第一项和最后一项减半. 2.3.2 龙贝格算法[10]在龙贝格算法中使用上述公式.设(,0)R n 表示具有2n个子区间的梯形估计,我们有[]1211(0,0)()()()21(,0)(1,0)((21))2n n n i R b a f a f b R n R n hf a i h -=⎧=-+⎪⎪⎨⎪=-++-⎪⎩∑ , (2.12) 对于一个适度的M 值,计算(0,0),(1,0),(2,0),,(,0)R R R R M L ,并且其中没有重复的函数值的计算.在龙贝格算法的其余部分中,还要计算附加值(,)R n m .所有这些都可以被理解为积分I 的估计.计算出(,0)R M 后,不再需要被积函数f 值的计算.根据公式[]1(,)(,1)(,1)(1,1)41m R n m R n m R n m R n m =-+-----, (2.13)对于1n ≥和1m ≥构造R 阵列的各列.定理 2.4(龙贝格算法收敛性定理)[10]若[],f C a b ∈,则龙贝格阵列中每一列都收敛于f 的积分.因此,对每个m ,lim (,)()baR n m f x dx =⎰.2.4高斯求积[11]前面研究的求积公式都是事先确定了n 个节点,然后按使求积公式阶数达到最大的原 则选取最佳权.由于自由参数为n 个,所以阶数一般为n-1,但如果节点的位置也自由选择,则自由参数的个数将变为2n ,因此求积公式的阶数可达到2n-1.高斯求积公式就是通过选择最佳的节点和权,使求积公式的阶数最大化.一般地,对每个n ,n 点高斯公式都是唯一的,而且阶数为2n-1.因而,对一定的节点个数,高斯求积公式的精度是最高的.但它的求得比牛顿—柯特斯公式要困难得多.虽然它的节点和权也可由待定系数法确定,但得到的方程是非线性的.2.4.1 高斯求积公式[11]为说明高斯求积公式,推导区间[]1,1-上的两点公式1112221()()()()()I f f x dx w f x w f x G f -=≈+=⎰,其中的节点1x 、2x 及权1w 、2w 按使求积公式阶数最大化的原则选取.令公式对前四个单项式精确成立,得力矩方程组112111122112221122113331122112,0,2,30.w w dx w x w x xdx w x w x x dx w x w x x dx ----⎧+==⎪⎪+==⎪⎪⎨⎪+==⎪⎪+==⎪⎩⎰⎰⎰⎰这个非线性方程组的一个解为12121,1,x x w w =-===另一个解可通过改变1x ,2x 的符号而得到.这样,两点高斯求积公式为2()(G f f f =-+,阶数为3.另外,高斯求积公式的节点也可以由正交多项式得到.若p 是n 次多项式,且满足()0,0,,1,bk ap x x dx k n ==-⎰L 则p 与[],a b 区间上所有次数小于n 的多项式正交,容易证明:1. p 的n 个零点都是实的、单的,且位于开区间(,)a b .2. 区间[],a b 上以p 的零点为节点的n 点插值型求积公式的阶数为2n-1,是唯一的n 点高斯公式.定义2.2[12] 如果1n +个节点的求积公式()()()nbk k ak x f x dx A f x ρ=≈∑⎰(2.14)的代数精度达到21n +,则称式(2.14)为高斯型求积公式,此时称节点k x 为高斯点,系数k A 称为高斯系数.定理2.5[12] 以01,,,n x x x L 为高斯点的插值型求积公式具有21n +次代数精确度的充要条件是以这些节点为零点的多项式101()()()()n n x x x x x x x ω+=---L与任意次数不超过n 的多项式()p x 带权()x ρ均在区间[],a b 上正交,即1()()()0bn ax p x x dx ρω+=⎰. (2.14)定理2.6 高斯公式()()nbi i ai f x dx A f x =≈∑⎰(2.15)的求积系数k A 全为正,且 2()(),0,1,,bbk k k aaA l x dx l x dx k n ===⎰⎰L . (2.16)定理2.7 对于高斯公式(2.14),其余项为 (22)211()()()()(22)!b n n a R f f x x dx n ξρω++=+⎰ , (2.17) 其中[]101,,()()()().n n a b x x x x x x x ξω+∈=---L2.4.2 高斯—勒让德(Gauss-Legendre )公式[13] 对于任意求积区间[],a b ,通过变换22a b b ax t +-=+,可化为区间[]1,1-,这时11()()222bab a a b b af x dx f t dt --+-=+⎰⎰. 因此,不失一般性,可取1,1a b =-=,考查区间[]1,1-上的高斯公式 11()()ni i i f x dx A f x -==∑⎰. (4.5)我们知道,勒让德多项式1211111()(1)2(1)!n n n n n d L x x n dx+++++⎡⎤=-⎣⎦+, (4.6) 是区间[]1,1-上的正交多项式,因此,1()n L x +的n+1个零点就是高斯公式(4.5)的n+1个节点.特别地,称1()n L x +的零点为高斯点,形如(4.5)的高斯公式称为高斯—勒让德公式.以上这些公式中的节点和求积系数可查表得到. 2.4.3 高斯—哈米特求积公式(Gauss-Hermite )[14] Gauss-Hermite 求积公式2()0()()nx n k k k ef x dx f x ω∞--∞=≈∑⎰, (4.7)其余项为(22)1(().2(22)!n n n n R f f n ξ+++=+ (4.8)2.4.4 高斯—切比雪夫(Gauss-Chebyshev )求积公式[15] 区间为[]1,1-,权函数()x ρ=Gauss 型求积公式,其节点k x 是Chebyshev多项式1()n T x +的零点,即21cos (0,1,,)2(1)k k x k n n π⎡⎤+==⎢⎥+⎣⎦L ,而(0,1,,)1k A k n n π==+L于是得到1021cos 12(1)nk k f n n ππ-=⎡⎤+≈⎢⎥++⎣⎦∑⎰(4.9) 称为Gauss-Chebyshev 求积公式,公式的余项为 (22)2(1)2()(),(1,1)2(22)!n n n R f f n πηη++=∈-+ , (4.10) 这种求积公式可用于计算奇异积分.2.5 递推型高斯求积[10]高斯求积公式不具有递推性:当节点个数一定时,如果自由选择所有的节点和权以达到最高的阶数,则节点个数不同的公式一般没有公共节点,这意味着与一组节点对应的积分值,在用另外一组节点计算积分值时不能被利用.Kronrod 求积公式避免了这种工作量的增加,这类公式是对称的,n 点高斯公式n G 与2n+1个点Kronrod 公式21n K +对应.21n K +节点的约束条件为:以n G 的节点作为21n K +的节点,按求积公式达到最高阶数的要求确定21n K +中剩下的n+1个节点及2n+1个权(其中包括n G 的节点的权).这样,求积公式的阶数可达到3n+1,而真正2n+1个点高斯公式应该是4n+1阶的,所以精度和效率是一对矛盾.使用两个节点个数不同的求积公式的主要原因是可以用它们的差估计积分近似值的误差.使用Gauss-Kronrod 公式对时,若以21n K +的值作为积分的近似值,则一半基于理论,一半基于经验,可以得到关于误差的保守估计: 1.521(200)n n G K +-.Gauss-Kronrod 公式不仅有效地提供了较高的精度,还给出了可靠误差估计,所以它被认为是最有效的求积公式之一,并且构成了主要软件库中求积程序的基础,特别地,公式715(,)G K 已被普遍使用.三、 总结部分因为一些定积分的求解比较复杂,所以数值积分的理论与方法一直是计算数学研究的基本课题.各种定积分的数值计算方法的出现和发展,加快和简化了求解定积分的效率和步骤.以上主要介绍了各种数值积分的方法——牛顿-科茨求积公式,复合求积公式,龙贝格积分法,高斯求积公式等.每种方法都有各自的优缺点,针对不同的积分函数采用不同的方法,所以在实际计算时,要做适当的采取.相信随着理论分析和研究的日益深入,求定积分的数值计算方法将更加简单和完善,为我们的计算带来前所未有的方便,在数学领域也将会更上一层楼.四、参考文献[1] 孙志忠,吴宏伟,袁慰平,闻震初.计算方法与实习(第4版)[M].南京:东南大学出版社,2009,(2): 128~129.[2]Micheal T .Heath . 张威,贺华,冷爱萍译.科学计算导论(第2版)[M].北京:清华大学出版社,2005,(10): 396~297.[3]李桂成.计算方法[M].北京:电子工业出版社,2005,(10):186.[4] 现代应用数学手册编委会. 现代应用数学手册——计算与数值分析卷[M]. 北京:清华大学出版社,2005,(1): 163~168.[5] 林成森. 数值计算方法(上)[M]. 北京:科学出版社,2004,(5): 220~221.[6]冯康.数值计算方法[M].北京:国防工业出版社,1978,(12): 45~47.[7]孙志忠,袁慰平,闻震初.数值分析(第2版)[M].南京:东南大学出版社,2002,(1): 191~194.[8] (美)柯蒂斯F .杰拉尔德 帕特里克O .惠特莱. 应用数值分析(第7版)[M].北京:机械工业出版社,2006,(8): 222~225.[9]夏爱生,胡宝安,孙利民,夏凌辉.复化Simpson 数值求积公式的外推算法[J].军事交通学院学报.2006,第8卷(第1期): 66~68.[10](美)David Kincaid, Ward Cheney .王国荣,俞耀明,徐兆亮译.数值分析(原书第三版)[M].北京:机械工业出版社,2005,(9): 400~403.[11]M.T.Heath. Scientific Computing:An Introductory Survey, Sscond Edition[M].清华大学出版社.英文影印版. 2001,(10): 351~355.[12]封建湖,车刚明,聂玉.数值分析原理[M].北京:科学出版社,2001,(9): 111~114.[13]杨大地,涂光裕.数值分析[M].重庆:重庆大学出版社,,2006,(9): 139~142.[14] 黄明游,刘播,徐涛.数值计算方法[M].北京:科学出版社,2005,(8):137~138.[15]Jeffery J.Leader. Numerical Analysis and Scientific Computation[M].英文影印本.北京:清华大学出版社,2005,(8): 342~349。
基础实验二 定积分数值计算
一、实验目的
学习定积分的数值计算方法,理解定积分的定义,掌握牛顿-莱布尼兹公式。
二、实验材料
2.1定积分的数值计算
计算定积分⎰b a
dx x f )(的近似值,可将积分区间n 等分而得矩形公式
n a b n a b i a f dx x f n i b a ---+≈∑⎰=])
1([)(1 或 n a
b n a b i a f dx x f n i b a --+≈∑⎰=][)(1 也可用梯形公式近似计算
n
a b b f a f n a b i a f dx x f n i b
a -++-+≈∑⎰-=]2)()()([)(11 如果要准确些,可用辛普森公式
n a
b b f a f a b i a f n a b i a f dx x f n i n i b a 6)]()()2)21((4)(2[)(111-++--++-+≈∑∑⎰=-= 对于⎰1
0sin xdx ,矩形公式、梯形公式、辛普森公式的Mathematica 程序为
a=0;b=1;k=10; f[x_]:=Sin[x];
d=N[Integrate[f[x],{x,a,b}],k];(计算精确值)
s1[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,0,m-1}],k];(取小区间左端点的矩形公式) s2[m_]:=N[Sum[f[a+(i+1/2)*(b-a)/m]*(b-a)/m,{i,0,m-1}],k]; (取小区间中点的矩形公式)
s3[m_]:=N[Sum[f[a+i*(b-a)/m]*(b-a)/m,{i,1,m}],k]; (取小区间右端点的矩形公式) s4[m_]:=N[Sum[(f[a+i*(b-a)/m]+f[a+(i+1)*(b-a)/m])/2*(b-a)/m,{i,0,m-1}],k]; (梯形公
式)
s5[m_]:=N[(b-a)/m/6*((f[a]+f[b])+2*Sum[f[a+i*(b-a)/m],{i,1,m-1}]
+4*Sum[f[a+(i-1/2)*(b-a)/m],{i,1,m}]),k];(辛普森公式)
r1[m_]:=d-s1[m];r2[m_]:=d-s2[m];r3[m_]:=d-s3[m];r4[m_]:=d-s4[m];r5[m_]:=d-s5[m];(误差)
t=Table[{s1[m],r1[m],s2[m],r2[m],s3[m],r3[m],s4[m],r4[m],s5[m],r5[m]},
{m,100,1000,100}]
利用以上程序计算⎰10dx 、⎰10xdx 、⎰102
dx x ,并对几个公式比较。
2.2可积条件
如果函数)(x f 在区间[]b a ,上连续,则)(x f 在区间[]b a ,上可积。
反之不然。
2.3牛顿-莱布尼兹公式
设函数)(x f 在[]b a ,上连续,而且)(x F 是)(x f 的一个原函数,则有牛顿-莱布尼兹公式⎰-=b a a F b F dx x f )()()(。
函数
⎩⎨⎧≠==0001)(x x x f 在[]2,1-不连续、不存在原函数,但在[]2,1-上可积;函数⎩⎨⎧>≤-=--00cos sin 22)(11x x x x x x x g 在[]2,1-不连续,但在[]2,1-上可积、存在原函数
⎩⎨⎧>≤=-00sin )(122
x x x x x x G 。
此外函数⎩⎨⎧∉∈=Q x Q x x D 01)(处处不连续、不存在原函数,在任意区间(长度大于0)上不可积。
求原函数并验证牛顿-莱布尼兹公式的Mathematica 程序
f[x_]:=Sin[x];
Integrate[f(x),x](求不定积分)
F[x_]:=%(定义原函数)
d=NIntegrate[f(x),{x,a,b}](求定积分)
df=F[b]-F[a] (计算原函数的增量)
r=d-df
三、实验准备
认真阅读实验目的与实验材料后要正确地解读实验,在此基础上制定实验计划(修改、补充或编写程序,提出实验思路,明确实验步骤),为上机实验做好准备。
四、实验思路提示
3.1定积分的定义
先对一个函数,例如x sin 在区间[0,1],在程序中改变m (例如10=m 、100、10000)并适当扩展有效数字(例如10=k 、20、50),运行程序计算定积分的近似值,分析误差。
再考虑其它函数。
最后对几个公式比较。
3.2牛顿-莱布尼兹公式
先对一个函数,例如x sin 在区间[0,1], 运行程序计算。
再考虑其它函数,例如指数函数、分段连续函数、)(x D 。
分析可积条件及牛顿-莱布尼兹公式成立的条件。