1、数值积分用matlab实现东北大学秦皇岛分校数值计算课程设计报告数值积分及Matlab实现学 院数学与统计学院专 业信息与计算科学学 号5133117姓 名楚文玉指导教师张建波 姜玉山成 绩教师评语：指导教师签字： 2015年07月14日1 绪论在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题1.1 课题的背景对于定积分在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式可以计算定积分的值，但在很多情况下的原函数不易求出或非常复杂被积函数的原函数很难用初等函数表达出来，例如等；有的函数的原函数存在，但其
2、表达式太复杂，计算量太大，有的甚至无法有解析表达式因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一随着计
3、算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义国内外众多学者在数值积分应用领域也提出了许多新方法在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab软件编程的实现，应用于实际生活中1.2 课题的主要内容框架1.2.1 数值积分各求积公式简介简介牛顿-柯特
4、斯求积公式及其辛普森求积公式，龙贝格求积公式，高斯求积公式的基本理论基础和方法1.2.2 求积公式的代码实现通过理解各种数值积分求积公式的原理方法，通过Matlab软件编程，实现以上求积公式1.2.3 应用举例通过简单举例，自建一个相对简单和复杂的函数，用上面编写的Matlab源程序来解决实际问题，体会数值积分和Matlab的优势2 牛顿柯特斯公式及Matlab实现2.1 牛顿柯特斯公式的基本原理方法设将积分区间a, b划分为n等分，步长为，选取等距节点构造出的差值型求积公式， （2.1）称为牛顿-柯特斯公式，式中称为柯特斯系数根据， （2.2）引进变量代换，则有 （2.3）当n = 2时，此
5、时柯特斯系数为，相应的求积公式就是辛普森求积公式： （2.4）2.2 牛顿柯特斯公式的Matlab实现functionC, g = NCotes(a, b, n, m)% a，b分别为积分的上下限；% n是子区间的个数； % m是被调用第几个被积函数； % 当n=1时计算梯形公式；当n=2时计算辛普森公式，以此类推；I = n; h = (b - a) / i; z = 0; for j = 0 : i x(j + 1) = a + j * h; s = 1; if j = 0 s = s; elsefor k=1 : j s =s * k; endendr = 1; if i - j = 0
6、 r = r; else for k = 1 : (I - j) r = r * k; end end if mod(I - j), 2) = 1 q = -(I * s * r); else q = i * s * r; endy = 1; for k = 0 : i if k = j y = y * (sym(t) - k); endend l = int(y, 0 , i); C(j + 1)= l / q; z = z + C(j + 1)*f1(m, x(j + 1); end g=(b - a)*z3 复合求积公式及Matlab实现3.1 复合梯形公式的基本原理将区间a, b划分成
7、n等分，分点，,在每个子区间（）上采用梯形公式得： （3.1）记 （3.2）称式(3.2)为复合梯形公式3.2 复合梯形公式的Matlab实现function s = trapr1(f, a, b, n)% f表示被积函数；% a，b表示积分上下限；% n是子区间的个数；h = (b - a) / n; s = 0; for k = 1 : (n - 1) x = a + h * k; s = s + feval(f, x); end format long s = h*(feval(f, a) + feval(f, b) / 2 + h * s;3.3 复合辛普森求积公式的基本原理将区间a,
8、b分等分，在每一个子区间上采用辛普森公式，若记，则得 （3.3）记（3.4）称式(3.4)为复合辛普森求积公式3.4 复合辛普森求积公式的Matlab实现function s = simpr1(f, a, b, n) % f表示被积函数；% a，b表示积分上下限；% n是子区间的个数；h = (b - a) / (2 * n); s1 = 0; s2 = 0; for k = 1 : n x = a + h * (2*k - 1); s1 = s1 + feval(f, x); end for k = 1 : (n - 1) x = a + h * 2 * k; s2 = s2 + feval
9、(f, x); end s = h*(feval(f, a) + feval(f, b) + 4 * s1 + 2 * s2) / 3;4 龙贝格求积公式及Matlab实现4.1 龙贝格算法的基本原理由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果作为积分的近似值，其误差近似值为 （4.1）可以设想，如果用这个误差作为的一种补偿，即将 （4.2）作为积分的近似值，可望提高其精确度直接根据复化求积公式，不难验证 （4.3）这说明，将区间对分前后两次复化梯形公式的值，按式 （4.4）作线性组合恰好等于复合辛普森公式的值，它比更接近于近似值 同样，根据 （4.5）用于作线性组合会得
10、到比更精确的值，且通过直接验证可得 （4.6）再由 （4.7）用与作线性组合，又可得到比更精确的值，通常记为，即 （4.8）此式（4.8）就称为龙贝格求积公式 上述用若干积分近似值推算出更为精确的积分近似值得方法，称为外推法我们将序列分别称为梯形序列，辛普森序列，柯特斯序列和龙贝格序列由龙贝格序列求积的算法称为龙贝格算法具体步骤为：第一步：算出和的值，根据公式 （4.9）求出;第二步：将区间a,b分半，算出的值，并根据(4.3)和(4.9)式计算和;第三步：再将区间分半，算出及的值，并根据(4.3)和(4.9)式计算和，再有公式(4.4)求出；第四步：再将区间分半，计算，, ,并根据公式(4.
11、5)计算第五步：再将区间分半，类似上述过程计算，, ，重复以上步骤即可得到， ,一直到龙贝格序列中前后两项的绝对值差不超过给定的误差险为止4.2 龙贝格算法的Matlab实现function R, quad, err, h=romber(f, a, b, n, delta) % f表示被积函数；% a，b表示积分上下限；% n是子区间的个数% delta是误差限M = 1; h = b - a; err = 1 J = 0; R = zeros(4, 4); R(1, 1) = h * (feval(f, a) + feval(f, b) / 2 while (err delta) & (J
12、n) ) | (J Ncotes(0, 1, 1, 2) 回车得到；ans = 0.9270输入Ncotes（0, 1, 2, 2） 回车得到：ans = 0.94616.1.2 复合梯形公式和复合辛普森求积公式的计算结果解：建立一个M文件定义一个f.m函数function y = f(x) if x = 0y = 1; else y = sin(x) / xend输入trapr1(f, 0, 1, 10) 回车得到：ans = 0.9458输入simpr1(f, 0, 1, 10)ans =0.94616.1.3 龙贝格求积公式的计算结果（取误差不超过）解：建立一个M文件定义一个f.m函数f
13、unction y = f(x) if x = 0y = 1; else y = sin(x) / xend输入romber（f, 0, 1, 10，0.5*10(- 10)） 回车得到：quad = 0.9461ans = 0.9207 0 0 0 0 0.9398 0.9461 0 0 0 0.9445 0.9461 0.9461 0 0 0.9457 0.9461 0.9461 0.9461 00.9460 0.9461 0.9461 0.9461 0.94616.1.4 高斯勒让德求积公式的计算结果(给定节点3)解：建立一个M文件定义一个f.m函数function y = f(x) ；
14、y = sin(x) / x;输入gauss8（f,0,1,-0.7745966692,0.7745966692,0,0.5555555556,0.5555555556,0.8888888888）回车得到：ans = 0.89566.2 复杂数值积分的解(精确值-1.8785)6.2.1 牛顿柯特斯当n = 1时的梯形算法和n = 2时的辛普森算法的结果解：先用M文件定义一个f1.m的函数function f = f1(i, x) g(1) = sqrt(x);if x = 0 g(2) = 1; else g(2) = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4
15、+ x.2);endf = g(i); 输入Ncotes(0, 2*pi, 1, 2)回车得到ans = -1.8692输入Ncotes(0, 2*pi, 2, 2) 回车得到ans = -1.87046.2.2 复合梯形公式和复合辛普森求积公式的计算结果解：先用M文件定义一个f1.m的函数function y = f(x) if x = 0y = 1; else y = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4 + x.2);end输入trapr1(f, 0, 2*pi, 10) 回车得到：ans = -1.8748输入simpr1(f, 0, 2*pi, 10
16、)ans = -1.86946.2.3 龙贝格求积公式的计算结果（取误差不超过）解：先用M文件定义一个f1.m的函数function y = f(x) if x = 0y = 1; else y = cos(x) 1 / (1 + x.2) 1 / (4.*sqrt(4 + x.2);end输入romber（f, 0, 2*pi ,10，0.5*10(-10)） 回车得到：quad = -1.87646.3 各个求积公式的比较分析（以的各个积分结果为例） 表6.1各个数值积分的比较积分方法牛顿柯特斯（梯形）牛顿-柯特斯（辛普森）复合梯形复合辛普森龙贝格高斯-勒让德精确值0.94610.9461
17、0.94610.94610.94610.9461实际值0.92700.94600.94580.94600.94600.8956误差0.1890.00010.00030.00010.00010.0505代数精度1233105牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级复合梯形积分方法比单独的梯形积分精度高，龙贝格方法收敛速度快、计算精度较高，但是计算量较大高斯求积方法积分精度高、数值稳
18、定、收敛速度较快，但是节点与系数的计算较麻烦，而且要求已知积分函数结 论本文主要讨论了数值积分的计算方法并通过MATLAB软件编程实现，通过前面的研究我们知道求数值积分近似值的计算方法很多，包括牛顿-柯特斯求积公式、复合求积公式、龙贝格求积公式、高斯求积公式等等其中牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级因此，一般情况下，代数精度越高，积分公式计算精度也越高但是高阶的牛顿-柯特斯
19、方法的收敛性没有保证，因此，在实际计算中很少使用高阶的牛顿-柯特斯公式复合梯形积分方法比单独的梯形积分精度高，它的积分精度和被积函数有关，还和复合积分时的步长有关复合辛普森积分公式比单独的辛普森积分公式高近7个数量级，效果明显龙贝格方法收敛速度快、计算精度较高，但是计算量较大高斯求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦、而且要求已知积分函数一般来说，牛顿-柯特斯方法的代数精度越高，数值积分的效果越好、越精确当积分区间比较大的时候，可以采用复化积分方法可以得到更好的效果；变步长积分方法不仅可以很好地控制计算误差，并且可以寻找到适当的积分步长；龙贝格积分方法可以更好地利用变步长复合积分公式得到的积分序列从而得到更为精确的数值结果，是一个较好的数值积分方法高斯求积方法精确度高，收敛性快也是一种很优秀的数值积分方法参考文献1 张德丰. Matlab数值分析与应用M. 北京: 国防工业出版社, 2007. 1-12 胡良剑, 孙晓君.MATLAB数学实验M. 北京: 高等教育出版社, 20063 李庆扬, 王能超, 易大义. 数值分析M. 北京: 清华大学出版社, 2008
