数值积分用matlab实现

东北大学秦皇岛分校数值计算课程设计报告数值积分及Matlab实现学院数学与统计学院专业信息与计算科学学号*******姓名楚文玉指导教师张建波姜玉山成绩教师评语：指导教师签字：2015年07月14日1 绪论在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题． 1.1 课题的背景对于定积分()ba f x dx ⎰在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式()()()baI f x dx F b F a ==-⎰可以计算定积分的值，但在很多情况下的原函数()f x 不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来，例如2sin (),x x f x e x-=等；有的函数()f x 的原函数()F x 存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外，许多实际问题中的被积函数()f x 往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab 软件编程的实现，应用于实际生活中． 1.2 课题的主要内容框架1.2.1 数值积分各求积公式简介简介牛顿-柯特斯求积公式及其辛普森求积公式，龙贝格求积公式，高斯求积公式的基本理论基础和方法． 1.2.2 求积公式的代码实现通过理解各种数值积分求积公式的原理方法，通过Matlab 软件编程，实现以上求积公式． 1.2.3 应用举例通过简单举例，自建一个相对简单和复杂的函数，用上面编写的Matlab 源程序来解决实际问题，体会数值积分和Matlab 的优势．2 牛顿－柯特斯公式及Matlab 实现2.1 牛顿－柯特斯公式的基本原理方法设将积分区间[a, b]划分为n 等分，步长为b ah n-=，选取等距节点k x a kh =+构造出的差值型求积公式()0()()nn n k k k I b a C f x ==-∑， （2.1）称为牛顿-柯特斯公式，式中()n k C 称为柯特斯系数．根据()bk a A l x dx =⎰， 0,1,2,,.k n = （2.2）引进变量代换x a th =+，则有()0000(1)()!()!n kn n n n n kj j j kj kh t j Cdt t j dt b a k j nk n k -==≠≠--==----∏∏⎰⎰ （2.3） 当n = 2时，此时柯特斯系数为(2)16C =，(2)146C =，(2)216C =，相应的求积公式就是辛普森求积公式：()4()()62b a a b S f a f f b -+⎡⎤=++⎢⎥⎣⎦（2.4）2.2 牛顿－柯特斯公式的Matlab实现function[C, 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 : ix(j + 1) = a + j * h;s = 1;if j == 0s = s;elsefor k=1 : js =s * k;endendr = 1;if i - j == 0r = r;elsefor k = 1 : (I - j)r = r * k;endendif mod((I - j), 2) == 1q = -(I * s * r); elseq = i * s * r; end y = 1; for k = 0 : i if k ~= jy = y * (sym('t') - k); end endl = int(y, 0 , i); C(j + 1)= l / q;z = z + C(j + 1)*f1(m, x(j + 1)); endg=(b - a)*z3 复合求积公式及Matlab 实现3.1 复合梯形公式的基本原理将区间[a , b ]划分成n 等分，分点k x a kh =+，,0,1,b ak n n-=,在每个子区间[1,k k x x +]（0,1,1k n =-）上采用梯形公式得： 11110()()[()()]()2k kn n bx k k n ax k k h I f x dx f x dx f x f x R f +--+=====++∑∑⎰⎰（3.1）记11100[()()][()2()()]22n n n k k k k k h hT f x f x f a f x f b --+===+=++∑∑ （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 longs = h*(feval('f', a) + feval('f', b)) / 2 + h * s; 3.3 复合辛普森求积公式的基本原理将区间[a,b]分等分，在每一个子区间1,k k x x +]上采用辛普森公式，若记1/212k k x x h +=+，则得1111/210()()[()4()()]()6k kn n bx k k k n ax k k h f x dx f x dx f x f x f x R f +--++====+++∑∑⎰⎰（3.3）记1111/211/2001[()4()()][()4()2()()]66n n n n k k k k k k k k h hS f x f x f x f a f x f x f b ---+++====++=+++∑∑∑（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 : 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;4 龙贝格求积公式及Matlab 实现4.1 龙贝格算法的基本原理由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果2N T 作为积分I 的近似值，其误差近似值为21()3N N T T - （4.1） 可以设想，如果用这个误差作为2N T 的一种补偿，即将22241()341NN N N N T T T T T -+-=- （4.2） 作为积分的近似值，可望提高其精确度．直接根据复化求积公式，不难验证22241()341N N N N N N T T S T T T -=+-=- （4.3）这说明，将区间对分前后两次复化梯形公式的值，按式222211()()341N N N N N N I T T T T T T ≈+-=+-- （4.4）作线性组合恰好等于复合辛普森公式的值N S ，它比2N T 更接近于近似值． 同样，根据2222211(S )(S )1541N N N N N N I S S S S ≈+-=+-- （4.5） 用2N S 于N S 作线性组合会得到比2N S 更精确的值，且通过直接验证可得222224S 1(S )1541N N N NN N S C S S -=+-=- （4.6） 再由222211(C )(C )6341N N N N N I C C C ≈+-=-- （4.7） 用2N C 与N C 作线性组合，又可得到比2N C 更精确的值，通常记为R N ，即322234C 1R (C )6341N N N NN N C C C -≈+-=- （4.8） 此式（4.8）就称为龙贝格求积公式．上述用若干积分近似值推算出更为精确的积分近似值得方法，称为外推法．我们将 序列{T },{S },{C },{R }N N N N 分别称为梯形序列，辛普森序列，柯特斯序列和龙贝格序列．由龙贝格序列求积的算法称为龙贝格算法．具体步骤为： 第一步：算出(a)f 和(b)f 的值，根据公式1211(a (2j 1))(N 2,k 1,2,)222N k NN j b a b a T T f N N-=--=++-==∑ （4.9） 求出1T ;第二步：将区间[a,b]分半，算出()2a bf +的值，并根据(4.3)和(4.9)式计算2T 和1S ; 第三步：再将区间分半，算出()4b a f a -+及(3)4b af a -+⨯的值，并根据(4.3)和(4.9)式计算4T 和2S ，再有公式(4.4)求出1C ；第四步：再将区间分半，计算8T ，4S , 2C ,并根据公式(4.5)计算1R ． 第五步：再将区间分半，类似上述过程计算16T ，8S , 4C ，2R ．重复以上步骤即可得到1R ，2R,一直到龙贝格序列中前后两项的绝对值差不超过给定的误差险为止．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 = 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)；5 高斯－勒让德求积公式及Matlab 实现5.1 高斯－勒让德求积公式的基本原理在高斯求积公式()()()nbk k ak f x x dx A f x ρ=≈∑⎰（5.1）中，若取权函数()x ρ=1，区间[-1,1]，则得公式1-1()()nk k k f x dx A f x =≈∑⎰（5.2）我们知道勒让德多项式是区间[-1,1]上的正交多项式，因此，勒让德多项式1()n P x +的零点就是求积公式(5.2)的高斯点，形如(5.1)式的高斯公式特别的称为高斯-勒让德求积公式．如下表5.1所示为高斯-勒让德求积公式的节点数和系数． 5.2 高斯－勒让德求积公式的Matlab 实现 function quad = gauss8(f ,a, b, x, A) N = length(x); T = zeros(1, N);T = (a + b) / 2 + ((b - a) / 2) * x;quad = ((b - a) / 2) * sum(A. * feval('f', T));表5.1 高斯-勒让德求积公式的节点数和系数6 各个求积公式的应用举例与比较分析6.1 简单数值积分10sin xdx x的解(精确值0.946083070367183)6.1.1 牛顿－柯特斯当n=1时的梯形算法和n=2时的辛普森算法的结果解：先用M文件定义一个f1.m的函数function f = f1(i, x)g(1) = sqrt(x);if x == 0g(2) = 1;elseg(2) = sin(x) / x;endf = g(i);输入>>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;elsey = sin(x) / xend输入>>trapr1(‘f’, 0, 1, 10)回车得到：ans =0.9458输入>>simpr1(‘f’, 0, 1, 10)ans =0.946110 ）6.1.3 龙贝格求积公式的计算结果（取误差不超过10解：建立一个M文件定义一个f.m函数function y = f(x)if x == 0y = 1;elsey = sin(x) / xend输入>>romber（‘f’, 0, 1, 10，0.5*10^(- 10)）回车得到：quad =0.9461ans =0.9207 0 0 0 00.9398 0.9461 0 0 00.9445 0.9461 0.9461 0 00.9457 0.9461 0.9461 0.9461 00.9460 0.9461 0.9461 0.9461 0.9461 6.1.4 高斯－勒让德求积公式的计算结果(给定节点3) 解：建立一个M 文件定义一个f.m 函数 function y = f(x) ； y = sin(x) / x; 输入 >>gauss8（‘f ’,0,1,[-0.7745966692,0.7745966692,0],[0.5555555556,0.5555555556,0.8888888888]） 回车得到： ans = 0.89566.2 复杂数值积分2201(cos 1x dx x π-+⎰的解(精确值-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; elseg(2) = cos(x) – 1 / (1 + x.^2) – 1 / (4.*sqrt(4 + x.^2)); end f = 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;elsey = 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)ans =-1.869410 ）6.2.3 龙贝格求积公式的计算结果（取误差不超过10解：先用M文件定义一个f1.m的函数function y = f(x)if x == 0y = 1;elsey = 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 各个求积公式的比较分析（以10sin xdx x的各个积分结果为例）表6.1各个数值积分的比较牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级．复合梯形积分方法比单独的梯形积分精度高，龙贝格方法收敛速度快、计算精度较高，但是计算量较大．高斯求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦，而且要求已知积分函数．结论本文主要讨论了数值积分的计算方法并通过MATLAB软件编程实现，通过前面的研究我们知道求数值积分近似值的计算方法很多，包括牛顿-柯特斯求积公式、复合求积公式、龙贝格求积公式、高斯求积公式等等．其中牛顿-柯特斯方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，辛普森方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；牛顿-柯特斯积分方法的误差比辛普森积分精度高两个数量级．因此，一般情况下，代数精度越高，积分公式计算精度也越高．但是高阶的牛顿-柯特斯方法的收敛性没有保证，因此，在实际计算中很少使用高阶的牛顿-柯特斯公式．复合梯形积分方法比单独的梯形积分精度高，它的积分精度和被积函数有关，还和复合积分时的步长有关．复合辛普森积分公式比单独的辛普森积分公式高近7个数量级，效果明显．龙贝格方法收敛速度快、计算精度较高，但是计算量较大．高斯求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦、而且要求已知积分函数．一般来说，牛顿-柯特斯方法的代数精度越高，数值积分的效果越好、越精确．当积分区间比较大的时候，可以采用复化积分方法可以得到更好的效果；变步长积分方法不仅可以很好地控制计算误差，并且可以寻找到适当的积分步长；龙贝格积分方法可以更好地利用变步长复合积分公式得到的积分序列从而得到更为精确的数值结果，是一个较好的数值积分方法．高斯求积方法精确度高，收敛性快也是一种很优秀的数值积分方法．参考文献[1]张德丰. Matlab数值分析与应用[M]. 北京: 国防工业出版社, 2007. 1-1[2]胡良剑, 孙晓君.MATLAB数学实验[M]. 北京: 高等教育出版社, 2006[3]李庆扬, 王能超, 易大义. 数值分析[M]. 北京: 清华大学出版社, 2008。