Matlab 数值积分 差值 插值 数据积分与数值微分.doc

第十四章 数据积分与数值微分 我们在高等数学里已完整地学习了求函数的积分和微分问题。但我们知道，并不是所有的积分与微分问题，都可以高等数学方法解决的。最常见的是函数列表函数时，求导与求微公式都不能使用，此时必然借助于数值积分与数值微分。这一章的内容与微分方程的求解有密切的联系。14.1 数值积分的基本概念14.1.1 求积公式要解决定积分的近似计算问题，首先要建立所谓的求积公式。 设是的函数，对于定积分：，如果被积函数是用解析式表示的连续函数，从理论上讲可以使用牛顿-莱布尼兹 (Newton-Leibniz) 公式，即： (14.1.1)但如果：(1) 找不到用初等函数表示的原函数，如等等。(2) 虽然找到了原函数，但因表达式过于复杂而不便于计算。因此，建立有效的数值积分的方法就很有必要了。从上所述，数值求积公式应该避免用原函数表示，最好能由被积函数的直接决定，这种做法是否可行呢？图14.1.1所示表示的是上与轴之间的面积。图14.1.1 积分示意图 由定积分的定义： (14.1.2)若等式的右边取某一近似值，即： (14.1.3)式中为求积节点，而称为求积系数，它仅与节点的选取有关，而与函数的具体形式无关，。上式便称为求积公式。 这类数值积分的方法通常称为机械求积法，它利用一些离散点上的函数的值作线性组合而得出积分的近似值。于是，求积分的问题转化为计算被积函数在节点处函数值的问题了。对于形如(14.1.3)的求积公式，关键在于确定求积节点和相应的求积系数14.1.2 求积公式的余项与代数精度求积公式与精确值的差为称求积公式的余项，也称截断误差为： (14.1.4)对于一个不超过m次的多项式，若其求积公式都能够精确等于定积分的数值，即，而对于m+1以上的多项式的却不能精确等于定积分的数值，则称该求积公式的代数精度为m。任何一个求积公式的代数精度至少为0，即取，求积公式精确地等于： (14.1.5)上式是求积系数应满足的基本条件，可以用它来检验一个求积公式的系数的正确性。一般地，要求(14.1.3)的代数精度为，只要令它对于都能精确成立，也就是说要求： (14.1.6)14.2 插值型求积公式14.2.1 插值型求积公式设给定一组节点及函数在这些节点上的值，按前面所讲的插值法，构造插值多项式。这样，其对应插值型求积公式为： (14.2.1)式中，求积系数通过插值基函数积分得出，即： (14.2.2)由插值余项定理，可得到插值型求积公式的余项为： (14.2.3)其中与变量有关，且。如果求积公式(14.2.1)是插值型的，由(14.2.3)式知，对于次数小于的多项式，其余项等于0，所以，此时求积公式的代数精度至少为。反之，如果的求积公式的代数精度至少是，则它必定是插值型的。14.2.2 内插求积公式由插值公式，对任一函数(包括离散形式的离散数据)，可用一个次多项式对其进行插值，当插值多项式为拉格朗日多项式时，即： (14.2.4)内插求积公式： (14.2.5)其中： (14.2.6) (14.2.7)14.2.3 牛顿-柯特斯公式牛顿-柯特斯 (Newton-Cotes) 公式：等距节点的内插求积公式。将积分区间等分，得牛顿牛顿-柯特斯公式： (14.2.8)其中： (14.2.9) 这种等距节点的内插求积公式通常称为牛顿-柯特斯公式。下面先介绍几个牛顿-柯特斯公式的特殊形式。14.2.4 梯形公式在牛顿-柯特斯公式中取得梯形公式为： (14.2.10)梯形公式是用梯形面积近似代替曲边梯形的面积。梯形公式的误差估计为：，() (14.2.11)梯形求积公式的代数精度为1。图14.2.1 梯形公式示意图14.2.5 抛物形公式在牛顿-柯特斯公式，取得抛物形公式 (又称辛普生Simpson公式)： (14.2.12)式中： ， (14.2.13)上式也可以写成： (14.2.14)抛物形公式的几何意义：抛物形公式是用抛物线围成的曲边梯形的面积近似代替围成的曲边梯形的面积。图14.2.2 抛物形公式示意图抛物形公式的误差估计为：，() (14.2.15)抛物形公式的代数精度为3，这是因为，若为三次以下多项式，式(14.4.14)精确成立。14.2.6 抛物形公式的MATLAB实现 程序功能：用抛物线公式逼近积分，其中。*function Z=srule(f,a0,b0)% f是被积函数。% a0和b0是被积区间的上下限。% 输出的值是所求的积分值。h=(b0-a0)/2;C=zeros(1,3);C=feval(f,a0),feval(f,(a0+b0)/2),feval(f,b0);S=h*(C(1)+4*C(2)+C(3)/3;Z=S;* 例14.2.1 用上述程序计算积分的值。解：先用m文件先定义一个名为f1421.m的函数文件。 function y = f1421(x) y = 1/x(1/2);建立一个主程序prog1421.m clcclear allfixpt(f1421,0.04,1) 然后在MATLAB命令窗口运行上述主程序，即： prog1421计算结果如下。ans = 1.8475这说明：。14.2.7 布尔公式当时，称为布尔公式：(14.2.16)其中，。14.2.8 牛顿-柯特斯公式牛顿-柯特斯公式的求积系数为： (14.2.17)其中，不依赖函数和积分区间，称为牛顿-柯特斯系数。表14.2.1 牛顿-柯特斯系数值123456牛顿-柯特斯公式的余项与代数精度：当，在上存在时，牛顿-柯特斯公式的余项为：，() (14.2.17)其中：。当是不超过次的多项式，牛顿-柯特斯公式的代数精度至少为，且当为偶数时，牛顿-柯特斯公式的代数精度可达到。牛顿-柯特斯公式 的几何意义：左图为上的梯形公式求积分。右图为上的抛物形公式求积分。 左图为上的辛普生公式求积分。右图为上的布尔公式求积分。14.3 复化求积公式 从求积系数的过程中可以看到，插值节点的增多 (的增大)，使用牛顿-柯特斯公式时，有可以导致求积系数出现负数，当时，牛顿-柯特斯求积系数会出现负数。另一方面，从求积公式的余项的讨论中看到，被积函数所用的插值多项式次数越高，对函数光滑性的要求也越高。如果我们将积分区间划分成若干个小区间，在各小区间上采用低次的求积公式 (梯形公式或抛物形公式)，然后再利用积分的区间可加性，把各区间上的积分加起来，得到新的求积公式。就是复化求积公式的思想。14.3.1 复化梯形公式将等分积分区间，分点为，其中。在每一个小区间是上用梯形公式，则有： (14.3.1) 如果将区间等分，这时节点数，除两端点外的节点记为，在每个小区间上再用梯形求积公式，则应为原节点和新增节点函数值的和的两倍加上两端点的函数值再乘上新区间长度的一半，即： (14.3.2)其中：。复化梯形公式的误差估计：设为二阶连续可微函数，则，为复化求积公式的余项： (14.3.3)其中。14.3.2 复化梯形公式的MATLAB程序程序功能：通过的个等步长节点逼近积分：其中，。*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);ends=h*(feval(f,a)+feval(f,b)/2+h*s;* 例14.3.1 用上述程序计算积分的值。解：先用m文件先定义一个名为f1431.m的函数文件。 function y = f1431(x) y = 1/(1+x/2);建立一个主程序prog1431.m clcclear allfixpt(f1431, -1, 1, 10) 然后在MATLAB命令窗口运行上述主程序，即： prog1431计算结果如下。ans = 1.5675即有：14.3.3 复化抛物形公式 复化抛物形求积公式为： (14.3.4)其中的下标表示将积分区间等分。复化抛物形公式的误差估计：设为连续四阶函数，则复化抛物形公式的余项： (14.3.5)其中。14.3.4 复化抛物形公式的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:n x=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;*例14.3.2 用上述程序计算积分的值。解：先用m文件先定义一个名为f1432.m的函数文件。 function y = f1432(x) y = 1/(1+x/2);建立一个主程序prog1431.m clcclear allfixpt(f1432, -1, 1, 10) 然后在MATLAB命令窗口运行上述主程序，即： prog1432计算结果如下。ans = 1.5708先定义一个函数文件，再调用上述函数运行。取，结果为：1.5708。14.4 龙贝格求积公式 数值积分公式计算的积分误差估计往往很复杂，涉及到被积函数的高阶微商。在处理实际问题时非常不便。因为函数的高阶微商根本无法求出 (如以数据表形式给出的函数)，或者十分复分。从前面可以看出，复化梯形公式比复化抛物形公式和复化Cotes公式精度低，收敛慢，这是其缺点，但是它的最大优点是算法简单，因此，我们关心的问题是如何发扬其优点，形成1个新的算法。这就是本节所讲的龙贝格 (Romberg) 求积法。采用龙贝格求积公式，逐次对分积分区间的方法，可以把前面计算的结果作为一个整体代入对分后的计算公式中，只需增加新的分点的函数值。龙贝格求积公式是一种很实用的公式。14.4.1 龙贝格求积公式怎样通过适当的线性组合，把复化梯形公式的近似值组合成更高精度的积分近似值呢？用复化梯形公式，取区间长度为，复化梯形公式的值与原积分值之间存在关系： (14.4.1)该式称为欧拉-麦克劳林求和公式。上式还表明，用近似的截断误差为： (14.4.2)龙贝格求积算法计算步骤如下：(1) 计算：(2) 对分区间并计算和：，其中为新分点的函数值之和。(3) 计算与：，(4) 利用外推公式：，()，()，直至求出。(5) 判断是否真，其中为给定的精度。若真，则，否则重复(2)、(3)、(4)的计算。表14.4.1 排成三角数表数表沿竖向和斜向都是收敛于，即当趋于无穷，与均收敛于。14.4.2 龙贝格求积公式的MATLAB程序程序功能：构造数表来逼近积分：，其中表示数表的最后一行，最后一列的值。*functionR,quad,err,h=romber(f,a,b,n,tol)% f是被积函数。% a, b分别为是积分的上下限。% n1是T数表的列数。% tol是允许误差。% R是T数表。% quad是所求积分值。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(errtol)&(Jn)|(J prog1441计算结果如下。quad = 3.1413ans =2.0000 0 0 0 0 0 0 02.7321 2.9761 0 0 0 0 0 02.9957 3.0836 3.0908 0 0 0 0 03.0898 3.1212 3.1237 3.1242 0 0 0 03.1233 3.1344 3.1353 3.1355 3.1355 0 0 03.1351 3.1391 3.1394 3.1394 3.1394 3.1394 0 03.1393 3.1407 3.1408 3.1408 3.1408 3.1408 3.1408 03.1408 3.1413 3.1413 3.1413 3.1413 3.1413 3.1413 3.1413由此可知：。14.5 高斯型求积公式14.5.1 最高代数精度的求积分式1. 最高代数精度简介高斯型求积公式讨论的就是最高代数精度的求积公式。下面讨论当节点数目固定为时，求积公式的最高代数精度能达到多少，又是通过怎样的方法达到最高代数精度。2. 最高代数精度的求法设求积公式的节点为：，求积公式为： (14.5.1)上式有个待定的求积数和个待定的节点，于是共有个待定系数。用，这时只需取，于是可解下列解方程组，求出求积系数： (14.5.2)式中：，() (14.5.3)所以，具有个节点的求积公式的代数精度最多为。另一方面，方程组(115.2)一定有解，也就是说求积公式(14.5.1)的代数精度可达到。3. 高斯型求积公式高斯型求积公式的求积系数：，()， (14.5.4)其中： (14.5.5)设，则高斯型求积公式的截断误差为： (14.5.6)其中：）14.5.2 几种常用的高斯型求积公式1. 高斯-勒让德（Gauss-Legendre）求积公式高斯-勒让德求积公式是古典的高斯求积公式，一般称为高斯求积公式。它是以积分区间上关于权函数的勒让德式多项式： (14.5.7)为正交多项式系，其中。表14.5.1列出了高斯-勒让德求积公式在时的节点与对应的系数。表14.5.1 高斯-勒让德求积公式的节点对应的系数20.5773503130.774596700.55555560.888888940.86113630.33998100.34785480.652145250.90617980.538469300.23692390.47862860.568888960.93246950.66120940.23861980.17132450.36076160.467913970.94910790.74153120.405845500.12948500.27970540.38183010.417959280.96028990.79666650.52553240.18343460.10122850.22238100.31370660.3626838高斯-勒让德求积公式的截断误差：， (14.5.8)2. 高斯-勒让德求积公式的MATLAB程序程序功能：利用高斯-勒让德求积公式求，使用变量替换：和，用高斯-勒让德求积公式的节点对应的系数表获得横坐标和系数，则有新节点，故有：。*function quad=gauss(f,a,b,X,A)% 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);*例14.5.1 用上述程序计算积分的值，给定节点为3。解：先用m文件先定义一个名为f1451.m的函数文件。 function y=f1451(x)y=sin(x)/x;建立一个主程序prog1451.m clcclearX=-0.7745966692, 0.7745966692, 0;A=0.5555555556, 0.5555555556, 0.8888888888;gauss(f1451, 0, 1, X, A) 然后在MATLAB命令窗口运行上述主程序，即： prog1451计算结果如下。ans = 0.8956即有： 3. 高斯-拉盖尔 (Gauss-Laguerre) 求积公式高斯-拉盖尔求积公式以积分区间上，关于权函数的拉盖尔多项式： (14.5.9)为正交多项式系，其中。高斯-拉盖尔求积系数为： (14.5.10)表14.5.2列出了高斯-拉盖尔求积公式在时高斯-拉盖尔求积系数。表14.5.2 相应节点xk的系数Ak20.58578643.41421360.85355340.146446630.41577462.29428046.28994510.71109300.27851770.010378940.32254771.74576114.53662039.39070910.60315410.35741870.03888790.000539350.26356031.41340313.59642587.085810012.64080080.52175560.39866680.07594240.00366180.0000234高斯-拉盖尔求积公式的截断误差为： (14.5.11)4. 高斯-埃尔米特 (Gauss-Hermite) 求积公式高斯-埃尔米特求积公式以积分区间上，关于权函数的埃尔米特多项式： (14.5.12)为正交多项式系。高斯-拉盖尔求积系数为： (14.5.13)表14.5.3列出了高斯-埃尔米特式多项式的零点与系数。表14.5.3 埃尔米特多项式的零点xk与系数Ak20.70710680.886226931.224744900.29540901.181635941.65068010.52464760.08131280.804914152.02018290