数值积分用matlab实现

东北大学秦皇岛分校数值计算课程设计报告数值积分及Matlab实现学院数学与统计学院专业信息与计算科学学号5133117姓名楚文玉指导教师张建波姜玉山成绩教师评语：指导教师签字：2015年07月14日数学与统计学院课程设计（实习）报告第1页1绪论在科研计算中，经常会碰到一些很难用公式定理直接求出精确解的积分问题，对于这类问题，我们一般转化为数值积分问题，用计算机来实现求解问题．1.1课题的背景对于定积分()bafxdx在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式()()()baIfxdxFbFa可以计算定积分的值，但在很多情况下的原函数()fx不易求出或非常复杂．被积函数的原函数很难用初等函数表达出来，例如2sin(),xxfxex等；有的函数()fx的原函数()Fx存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式．因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的．另外，许多实际问题中的被积函数()fx往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值．因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算．而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值．微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节数值积分是数学上重要的课题之一，是数值分析中重要的内容之一．随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域．现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义．国内外众多学者在数值积分应用领域也提出了许多新方法．在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等．通过这个课题的研究，我们将会更好地掌握运用数值积分算法求出特殊积分函数的定积分的一些基本方法、理论基础；并且通过Matlab软件编程的实现，应用于实际生活中．1.2课题的主要内容框架1.2.1数值积分各求积公式简介数学与统计学院课程设计（实习）报告第2页简介牛顿-柯特斯求积公式及其辛普森求积公式，龙贝格求积公式，高斯求积公式的基本理论基础和方法．1.2.2求积公式的代码实现通过理解各种数值积分求积公式的原理方法，通过Matlab软件编程，实现以上求积公式．1.2.3应用举例通过简单举例，自建一个相对简单和复杂的函数，用上面编写的Matlab源程序来解决实际问题，体会数值积分和Matlab的优势．2牛顿－柯特斯公式及Matlab实现2.1牛顿－柯特斯公式的基本原理方法设将积分区间[a,b]划分为n等分，步长为bahn，选取等距节点kxakh构造出的差值型求积公式()0()()nnnkkkIbaCfx，（2.1）称为牛顿-柯特斯公式，式中()nkC称为柯特斯系数．根据()bkaAlxdx，0,1,2,,.kn（2.2）引进变量代换xath，则有()0000(1)()!()!nknnnnnkjjjkjkhtjCdttjdtbakjnknk（2.3）当n=2时，此时柯特斯系数为(2)016C，(2)146C，(2)216C，相应的求积公式就是辛普森求积公式：()4()()62baabSfaffb（2.4）2.2牛顿－柯特斯公式的Matlab实现function[C,g]=NCotes(a,b,n,m)数学与统计学院课程设计（实习）报告第3页%a，b分别为积分的上下限；%n是子区间的个数；%m是被调用第几个被积函数；%当n=1时计算梯形公式；当n=2时计算辛普森公式，以此类推；I=n;h=(b-a)/i;z=0;forj=0:ix(j+1)=a+j*h;s=1;ifj==0s=s;elsefork=1:js=s*k;endendr=1;ifi-j==0r=r;elsefork=1:(I-j)r=r*k;endendifmod((I-j),2)==1q=-(I*s*r);elseq=i*s*r;endy=1;fork=0:i数学与统计学院课程设计（实习）报告第4页ifk~=jy=y*(sym('t')-k);endendl=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等分，分点kxakh，,0,1,baknn,在每个子区间[1,kkxx]（0,1,1kn）上采用梯形公式得：111100()()[()()]()2kknnbxkknaxkkhIfxdxfxdxfxfxRf（3.1）记11100[()()][()2()()]22nnnkkkkkhhTfxfxfafxfb（3.2）称式(3.2)为复合梯形公式．3.2复合梯形公式的Matlab实现functions=trapr1(f,a,b,n)%f表示被积函数；%a，b表示积分上下限；%n是子区间的个数；h=(b-a)/n;s=0;fork=1:(n-1)x=a+h*k;s=s+feval('f',x);数学与统计学院课程设计（实习）报告第5页endformatlongs=h*(feval('f',a)+feval('f',b))/2+h*s;3.3复合辛普森求积公式的基本原理将区间[a,b]分等分，在每一个子区间1,kkxx]上采用辛普森公式，若记1/212kkxxh，则得1111/2100()()[()4()()]()6kknnbxkkknaxkkhfxdxfxdxfxfxfxRf（3.3）记1111/211/2001[()4()()][()4()2()()]66nnnnkkkkkkkkhhSfxfxfxfafxfxfb（3.4）称式(3.4)为复合辛普森求积公式．3.4复合辛普森求积公式的Matlab实现functions=simpr1(f,a,b,n)%f表示被积函数；%a，b表示积分上下限；%n是子区间的个数；h=(b-a)/(2*n);s1=0;s2=0;fork=1:nx=a+h*(2*k-1);s1=s1+feval('f',x);endfork=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;数学与统计学院课程设计（实习）报告第6页4龙贝格求积公式及Matlab实现4.1龙贝格算法的基本原理由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果2NT作为积分I的近似值，其误差近似值为21()3NNTT（4.1）可以设想，如果用这个误差作为2NT的一种补偿，即将22241()341NNNNNTTTTT（4.2）作为积分的近似值，可望提高其精确度．直接根据复化求积公式，不难验证22241()341NNNNNNTTSTTT（4.3）这说明，将区间对分前后两次复化梯形公式的值，按式222211()()341NNNNNNITTTTTT（4.4）作线性组合恰好等于复合辛普森公式的值NS，它比2NT更接近于近似值．同样，根据2222211(S)(S)1541NNNNNNISSSS（4.5）用2NS于NS作线性组合会得到比2NS更精确的值，且通过直接验证可得222224S1(S)1541NNNNNNSCSS（4.6）再由222211(C)(C)6341NNNNNICCC（4.7）用2NC与NC作线性组合，又可得到比2NC更精确的值，通常记为RN，即322234C1R(C)6341NNNNNNCCC（4.8）此式（4.8）就称为龙贝格求积公式．上述用若干积分近似值推算出更为精确的积分近似值得方法，称为外推法．我们将序列{T},{S},{C},{R}NNNN分别称为梯形序列，辛普森序列，柯特斯序列和龙贝格序数学与统计学院课程设计（实习）报告第7页列．由龙贝格序列求积的算法称为龙贝格算法．具体步骤为：第一步：算出(a)f和(b)f的值，根据公式1211(a(2j1))(N2,k1,2,)222NkNNjbabaTTfNN（4.9）求出1T;第二步：将区间[a,b]分半，算出()2abf的值，并根据(4.3)和(4.9)式计算2T和1S;第三步：再将区间分半，算出()4bafa及(3)4bafa的值，并根据(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((errdelta)&(Jn))|(J4)J=J+1;h=h/2;s=0;数学与统计学院课程设计（实习）报告第8页forp=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;forK=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高斯－勒让德求积公式的基本原理在高斯求积公式0()()()nbkkakfxxdxAfx（5.1）中，若取权函数()x=1，区间[-1,1]，则得公式1-10()()nkkkfxdxAfx（5.2）我们知道勒让德多项式是区间[-1,1]上的正交多项式，因此，勒让德多项式1()nPx的零点就是求积公式(5.2)的高斯点，形如(5.1)式的高斯公式特别的称为高斯-勒让德求积公式．如下表5.1所示为高斯-勒让德求积公式的节点数和系数．5.2高斯－勒让德求积公式的Matlab实现functionquad=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));数学与统计学院课程设计（实习）报告第9页表5.1高斯-勒让德求积公式的节点数和系数nkxkAnkxkA00.00000002.000000030.86113630.33998100.34785480.652145210.57735031.000000040.90619780.53846930.00000000.23692690.47862870.568888920.77459670.00000000.55555560.888888950.93246950.66120940.23861920.17132450.36076160.46791396各个求积公式的应用举例与比较分析6.1简单数值积分10sinxd