第5章 数值积分.ppt

1、2020/9/3,1,数值计算与MATLAB,2,第五章 数值积分,我们知道,若函数f(x)在区间a,b上连续且其原函数为F(x),则可用Newton-Leibnitz公式 求定积分的值 , Newton- Leibnitz公式 无论在理论上还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：,3,(1) 被积函数f(x)并不一定能够找到用初等函数的 有限形式表示的原函数F(x)，例如： Newton-Leibnitz公式就无能为力了,(2) 还有被积函数f(x)的原函数能用初等函数表示， 但

2、表达式太复杂，例如函数,并不复杂， 但积分后其表达式却很复杂, 积分后其原函数F(x)为：,4,(3) 被积函数f(x)没有具体的解析表达式, 其函数 关系由表格或图形表示。 对于这些情况, 要计算积分的准确值都是十分困难的。由此可见, 通过原函数来计算积分有它的局限性, 因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题, 这时需要用数值解法来建立积分的近似计算方法。 将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想。,5,5.1计算积分的matlab符号法 使用matlab的符号计算功能，可以计算出许多积分的解析解和

3、精确解，只是有些精确解显示冗长繁杂，这时可以用vpa或eval函数把它转换成位数有限的数字，有效数字的长度可按需选取。 求积分的符号运算命令为int（取自integerate的前3个字母），调用格式为 s=int(fun,v,a,b) 1)输入参量fun是被积函数的符号表达式，可以是函数向量或函数矩阵。 2）输入参量v是积分变量，必须被界定成符号变量。如果被积函数中只有一个变量时可以缺省。,6,3）输入参量a、b为定积分的积分限，缺省时输出被积函数fun的一个原函数。 4）输入参量s为积分结果。s为有理表达式有时过于冗长时，可用vpa或eval命令把它转换成有限长度的小数。 例5.1 计算,7

4、,解：该函数图线的生成方法如下： y=1:10; ff=exp(-y.2)+log(y); plot(ff,linewidth,2),title(exp(-y.2)+log(y),grid,8,在matlab命令窗口中输入 int(exp(-y2)+log(y),1,10) 回车得到 ans = 1/2*pi(1/2)*erf(10)-9+10*log(2)+10*log(5)-1/2*pi(1/2)*erf(1) 输出的结果比较复杂，下面是用两种方式进行转换的输出结果 是比较他们的差别 输入 eval(1/2*pi(1/2)*erf(10)+10*log(2)+10*log(5)-9-1/2

5、*pi(1/2)*erf(1) ans = 14.1653,9,输入 vpa(1/2*pi(1/2)*erf(10)+10*log(2)+10*log(5)-9-1/2*pi(1/2)*erf(1) ans = 14.165253722580787828429530852382 例5.2 计算,10,解：该函数的图像生成方法如下 x=-10:1:10; ff=1./(x.2+1); plot(ff,linewidth,2),title(1./(x.2+1),grid,11,在matlab命令窗口中输入 syms x; f=1/(x2+1); a=int(f,-inf,inf) 回车得到 a =

6、 pi 得到精确值,12,例5.3计算 解：在matlab命令窗口中输入 syms x t; f=2*t*x; a=int(f,x,1,sin(t) 回车得到 a = t*(sin(t)2-1) 例5.4求 ，准确结果是,13,解：在matlab命令窗口中输入 syms x; Is=int(sqrt(log(1/x),x,0,1) 回车的到 Is = int(log(1/x)(1/2),x = 0 . 1) 输入 vpa(Is) ans = .88622692545275801364908374167057 再输入 vpa(Is,5) ans = .88623,14,例5.5计算单位域上的积分

7、 先把二重积分 转换为二次积分的形式 解在matlab命令窗口中输入 syms x y; Q=int(int(exp(-x2/2)*sin(x2+y),x,-sqrt(1-y2),sqrt(1-y2),y,-1,1) 回车得到 Q = int(-i*pi(1/2)*(exp(i*y)2*erf(1/2*(1-y2)(1/2)*(2-4*i)(1/2)*(2+4*i)(1/2)-erf(1/2*(1-y2)(1/2)*(2+4*i)(1/2)*(2-4*i)(1/2)/(2-4*i)(1/2)/exp(i*y)/(2+4*i)(1/2),y = -1 . 1) 在输入 vpa(Q,6) 回车得

8、ans = .536860-.171355e-8*i,15,例5.6计算 解：在命令窗口输入 syms x y; f=x2+y2; xlower=sqrt(y); xupper=2; ss=int(int(f,x,xlower,xupper),y,1,4); vpa(ss,6) 回车得到 ans = 9.58095,16,例5.7p176页 解 在matlab命令窗口中输入 syms x t,A=cos(x*t) exp(2*x*t) sin(x*t); -sin(x*t) 5*x2 -cos(x*t); int(A,t,-1,2) 回车得到 ans = (sin(2*x)+sin(x)/x,

9、 1/2*(exp(4*x)-exp(-2*x)/x, -(cos(2*x)-cos(x)/x (cos(2*x)-cos(x)/x, 15*x2, -(sin(2*x)+sin(x)/x,17,5.2 机械求积,1、数值求积的基本思想,由x=a,x=b,y=0,y=f(x)所围成的曲边梯形的面积,依据积分中值定理：,即：底为b-a，高为f()的矩形面积恰为积分值I,问题： 确定的位置， 求出f(）,18,数值求积方法：,对平均高度f()提供一种数值算法,若f()近似的用积分区间端点处的函数值f(a)和f(b)的算术平均值代替，便可导出计算积分的梯形公式：,若f()近似的用积分区间中点(a+b

10、)/2处的函数值代替，便可导出计算积分的中矩形公式：,19,若f()近似的用a,b,(a+b)/2三点高度的加权平均值1/6f(a)+4f(c)+f(b) ,c=(a+b)/2代替，便可导出计算积分的辛普森公式：,一般情况下的数值求积方法：,在积分区间a,b上适当选取某些节点xk,然后用f(xk)加权平均，得到f()近似值。 这样构造出来的求积公式为：,xk为求积节点 Ak为求积系数，亦称伴随节点xk的权,可以看出：求积系数Ak仅仅与节点xk的选取有关，而不依赖于被积函数f(x）的具体形式,20,2.代数精度的概念,研究近似方法结果的精度,对于一切次数m的多项式是准确的，但对于m+1次多项式不

11、准确，则称它具有m次代数精度,一般地，欲使该求积公式具有m次代数精度，只要令它对于f(x)=1,x,x2,xm都准确成立即可,21,梯形公式：,所以，梯形公式具有一次代数精度,22,中矩形公式：,所以，中矩形公式具有一次代数精度,23,辛普森公式：,所以，辛普森公式具有三次代数精度,24,25,可用代数精度作为标准来构造求积公式：,26,一般情况：,对于给定的一组求积节点xk(k=0,a,n),可以确定相应的求积系数Ak,使求积公式至少具有n次代数精度,求积公式的构造本质上是个解线性方程组的代数问题,27,3.插值型的求积公式,设已给出f(x)在节点xk(k=0,1,n)的函数值，作插值多项式

12、,由于pn(x)的求积是容易的，,此求积公式称为插值型的,28,对于任意次数n的多项式f(x),其插值多项式pn(x）就是它自身，因此此插值型的求积公式至少有n次代数精度,反之,如果该求积公式至少具有n次代数精度，则它对于插值基函数lk(x）是准确成立的,定理1：形如 的求积公式至少具有n次代数精度的充分必要条件是，它是插值型的。,这样，在求积节点xk已给出的情况下，求积系数Ak的确定就有两种可选择的方法 解线性方程组 计算积分,29,例1,30,(2),31,5.3 Newton-Cotes求积公式,截断误差为,32,33,表6-1,34,Newton-Cotes求积公式的代数精度,35,3

13、6,梯形公式：,37,Simpson公式：,38,39,如果积分区间比较大，直接地使用上述求积公式， 精度难以保证。,高次插值有Runge 现象，故采用分段低次插值 分段低次合成的 Newton-Cotes 复合求积公式。,（1）等分求积区间，比如取步长 ，分a, b为n等分，,分点为 k = 0, 1, 2, n,（2）在区间 xk, xk+1上使用以上求积公式求得Ik,（3）取和值 ，作为整个区间上的积分近似值。,5.5 复合求积公式,40, 复合梯形公式：,在每个 上用梯形公式：,= Tn,/*积分中值定理*/,41,复合梯形求积公式的matlab实现 function s=trapr1

14、(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); end s=h*(feval(f,a)+feval(f,b)/2+h*s;,42,例 利用程序计算积分 解：先用M文件定义一个名为f.m的函数 function y=f(x) y=1/(1+x2); 在matlab命令窗口中输入 trapr1(f,-1,1,10) 回车得到 ans = 1.5675,43, 复化 Simpson 公式：,= Sn,44,复合抛物形求积公

15、式的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); end for k=1:(n-1) x=a+h*2*k; s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;,45,例 利用程序计算积分 解：先用M文件定义一个名为f.m的函数 function

16、y=f(x) y=1/(1+x2); 在matlab命令窗口中输入 simpr1(f,-1,1,10) 回车得到 ans = 1.5708,46,一、梯形法的递推化逐次分半法,上一节介绍的复合求积方法对提高精度是行之有效的，但在使用求积公式之前必须给出合适的步长，步长取得太大精度难以保证，步长太小则会导致计算量的增加，而事先给出一个恰当的步长又往往是困难的 实际计算中常常采用变步长的计算方案，即在步长逐次分半(即步长二分)的过程中，反复利用复合求积公式进行计算，直至所求得的积分值满足精度要求为止,设将求积区间a，b分成n等分，则一共有n+1个分点，按梯形公式计算积分值Tn，需要提供n+1个函数

17、值如果将求积区间再二分一次，则分点增至2n+1个，我们来考察二分前后两个积分值之间的联系,5.7 龙贝格求积公式,逐次分半计算方案的实现:,47,注意到每个子区间xk，xk+1经过二分只增加了一个分点 xk+1/2( xk+xk+1)/2，用复合梯形公式求得该子区间上的积分值为,这里 代表二分前的步长.将每个子区间上的积分值相加得,48,二、龙贝格算法,根据复化梯形公式的余项表达式,可见，利用两种步长计算的结果能估计截断误差.若将该截断误差加到计算结果中，,就得出“改进的梯形求积公式”：,事后误差估计,49,这就是说用梯形法二分前后的两个积分值Tn与T2n的线性组合的结果得到复化辛普森法求积公

18、式,类似的情况，用辛普森法二分前后的两个积分值Sn与S2n的线性组合的结果可得到复化柯特斯求积公式,重复同样的手续，用柯特斯法二分前后的两个积分值Cn与C2n的线性组合的结果可得到龙贝格(Romberg)求积公式,我们在变步长的过程中运用加速公式，就能将粗糙的梯形值Tn逐步加工成精度较高的辛普森值Sn 、柯特斯值Cn和龙贝格值Rn .,一般有：, Romberg 算法：, ?, ?, ?, ,Romberg 序列,50,练习：用龙贝格求积法计算积分 的近似值.,51,龙贝格求积公式的matlab实现,构造T数表来逼近积分 其中， 表示T数表的最后一行，最后一列的值。,52,function R

19、,quad,err,h=romber(f,a,b,n,delta) % f是被积函数; % a,b分别是积分的上下限; % n+1是T数表的列数 % delta是允许误差 % 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 (errdelta) end quad=R(J+1,J+1),53,例 用程序计算 的值。精确到小数点后第10位。 romber(sqrt(4*x-x2,0,2,7,0.5*(10(-10) 回车得 quad = 1.1071 a

20、ns = 1.2000 0 0 0 0 0 0 0 1.1000 1.0667 0 0 0 0 0 0 1.1038 1.1051 1.1077 0 0 0 0 0 1.1063 1.1071 1.1073 1.1073 0 0 0 0 1.1069 1.1071 1.1071 1.1071 1.1071 0 0 0 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071 0 0 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071 0 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071 1.1071,54,5.8高斯公式,（1）含义： 积分公式的一般形式； 以前的节点是按等间距来选择，为了获得更高的代数精度节点也可以作为待定值。 （2）一点高斯公式 设一点高斯公式的形式为： 其实 都是需要待定的值。根据代数精度概念， 令 ，使积分公式准确成立，有,55,解得： ， 故一点高斯公式为： ，即为中矩形公式，它具有1次代数精度。 （3）二点高斯公式 设二点高斯公式的形式为： 其实 都是需要待定的值。根据代数精度概念， 令 ，使积分公式准确成立，有,56,该方程组不是线性方程组，故其求解比较困难