Chapt-5 数值微积分的数值解法.ppt

1、第5章 微积分的数值解法,在现代工程领域或者科研过程中所遇到的许多实际问题的求解，常常归结为定积分的计算。例如，力学和电气中功和功率的计算、电流电压平均值及有效值的计算，几何学中面积和何种及其重心的计算，化学化工中的热容、热焓、转化率、反应时间和反应器体积等等的计算，都与某些定积分的计算有关。但很多情况下，其积分函数f(x)的关系难以明确表示，即使能用解析明确表示，有时表达式也很复杂，不能用于实际计算。本章介绍计算积分和导数的实用数值方法。,5.1 数值积分的基本思想,一、定积分的数学含义 由定积分的定义： 其中 为第i个积分区间 内任一点。当n足够大时有以下近似表达式： 将被积函数在积分区间

2、上足够多的离散点的函数值与其所在小区间长度乘积之和作为定积分的近似数值。,(5-1),(5-2),5.1 数值积分的基本思想,二、定积分的几何含义 （5-2）式告诉我们，定积分的值就是由下面四条曲线所转面的面积，即：,图5-1 定积分的几何意义,5.1 数值积分的基本思想,如图5.1中所示的阴影部分的面积。这就是定积分的几何意义。,当求积区间只有一个时 1. 由积分中值定理，在（a, b）中必存在一点，使得 2.用f(a)和f(b)的平均值近似代替f()，则： 3.用中点的函数值近似代替f()，则：,5.1 数值积分的基本思想,推广到更一般情况，当把求积区间(a,b)分成n个小求积区间时，则根

3、据定积分的几何定义，我们有 其中 为小区间长度 分点为 ，它们所对应的纵坐标（函数）分别为 。如图5-2，取其中任一个小区间 进行分析。该区间为一曲边梯形，y用以下方法求其面积。,5.1 数值积分的基本思想,1. 用小矩形或小梯形近似表示时，然后求和得到I , 这就是高等数学中的矩形积分法和梯形积分法。它们都是用直线（一次函数）代替曲线（f(x)函数）求得的近似值，精度较低，一般要求N取得足够大才能保证精度。 2. 当用曲线去逼近f(x)时，只要曲线近似程度足够高，一般都能满足计算精度要求。 本章按照上述思想，主要介绍梯形积分法、抛物线积分法（Simpson法)、Newton-Cotes积分法

4、，龙贝格积分法以及高斯积分法,5.2 梯形积分法,一 定步长梯形积分法 1. 算法分析 设在区间(a, b)上有可积函数f(x),求积分值 将区间(a, b)分成n个相等的小区间，每个小区间之长为： h=(b-a)/n 各分点： x0=a, x1, x2, , xi, xi+1, , xn=b 其中 xi=a+i*h (i=1,2,3, ,n) 用p(x)=c*x+d直线近似代替f(x)(参见图5-2）,5.2 梯形积分法,用插值的方法，我们可求得 将其代入积分公式有,5.2 梯形积分法,第i个区间： Ti=(f(xi)+f(xi+1)*h/2 第i-1个区间： Ti-1=(f(xi-1)+f

5、(xi)*h/2 求和得到： 定步长梯形积分法的截断误差为:,(5-3),5.2 梯形积分法,该误差按h的2次方的速度下降,即定步长梯形积分法的误差阶为2，也就是说，定步长梯形积分法具有一阶代数精度（因为代数精度等于误差阶减1，所以2-1=1）。 2. 定步长梯形积分法的程序框图与通用程序设计 定步长梯形积分法的通用程序框图如图5-3所示，共分为三个程序框图，即在图中（1）为调用子程序的主控程序框图，（2）为函数子程序框图，（3）为积分通用子程序框图。,5.2 梯形积分法,（1）主控程序框图,（2）函数子程序框图,START,INPUT a,b,n,CALL sub T,OUTPUT sol,

6、END,FUNCTION F,F=expr,END FUNCTION,图5-3 梯形积分法的程序框图,5.2 梯形积分法,subroutine,h=(b-a)/n,Comp f(a),f(b),T=(f(a)+f(b)/2,DO i=1,n-1,x=a+i*h,Comp f(x),T=T+f(x),END DO i,T=T*h,END subroutine,（3）梯形积分子程序框图,5.2 梯形积分法,10C 主控程序 20 PROGRAM main 30 read(5,*) a,b,n 40 call SUBRPUTINE 50 T2=T1/2,Comp f(x),T2=T2+f(x),EN

7、D DO i,n=2*n;T1=T2,END subroutine,图5-4 梯形积分法子程序框图,n=1,x=a+(2*i-1)*h,abs(T2-T1)E?,Y,N,5.2 梯形积分法,10 SUBRPUTINE btxjf(a,b,n,T) 20 h=(b-a) 30 T1=h*(f(a)+f(b)/2 40 n=1 45 h=h/2;T2=T1/2 50 DO i=1,n 60 x=a+(2*i-1)*h 70 T2=T2+f(x)*h 80 END DO 90 if (abs(T2-T1)e) then 100 n=2*n;T1=T2 110 goto 45,120 endif 13

8、0 END SUBRPUTINE btxjf,5.2 梯形积分法,3.应用实例分析 例5.1 已知圆周率的近似值可以由下列定积分求出 试求出圆周率的似近值。 解： 因为 所以求解此问题的函数程序程序为 80C 函数子程序 90 FUNCTION f(x) 100 f=4.0/(1+x*x) 110 END FUNCTION f,5.2 梯形积分法,表5.1 梯形积分法输入不同N地的计算机输出结果,由上表可以看出，定步长梯形积分法的收敛速度是很慢的，对于本例，N=300时，即把积分区间0，1划分为300个小区间时才使计算结果精确到 0.00001。,5.3 Simpson积分法,梯形积分法是用直

9、线段（一次二项式）去逼近曲线f(x),而Simpson 是用抛物线段（二次三项式）去逼近曲线f(x). 一、定步长Simpson积分法 1.算法分析 设 将（a, b）分为n个（偶数）相等的小区间，则 求积节点坐标 对应求积点函数值,5.3 Simpson积分法,被积曲线用 y=f(x) 表达，而逼近曲线则采用 y=Ax2+Bx+C,5.3 Simpson积分法,第一段：x0,x1,x2 有 同理，第二段：x2,x3,x4 有 第末段： xn-2,xn-1,xn 有,求出系数A, B, C,5.3 Simpson积分法,共有n/2小块面积，对它们求和即可得Simpson法求积公式 或 这里,(

10、5-5),(5-6),5.3 Simpson积分法,同理我们可以推出辛普森积分法的截断误差为： 由上式可知辛普森积分法的误差阶和代数精度分别为4和3，在积分小区间总数与梯形积分法相同，辛普森积分法的精度较高。 2.通用程序框图与通用程序设计 这里只画出Simpson定步长积分法通用子程序框图（3）。,(5-7),5.3 Simpson积分法,subroutine,h=(b-a)/n,Comp f(a),f(b),S=f(a)-f(b),DO i=1,m,Comp f(x1), f(x2),S=S+4*f(x1)+ 2*f(x2),END DO i,x2=a+(2*i)*h,END subrou

11、tine,图5-5 Simpson法子程序框图,m=n/2,x1=a+(2*i-1)*h,S=S*h/3,5.3 Simpson积分法,10 SUBRPUTINE Simpson(a,b,n) 20 h=(b-a)/n 30 m=n/2 40 S=f(a)-f(b) 50 DO i=1,m 60 x1=a+(2*i-1)*h 70 x2=a+2*i*h 80 S=S+4.0*f(x1)+2.0*f(x2) 90 END DO 100 S=S*h/3.0 110 END SUBROUTINE Simpson,5.3 Simpson积分法,二、变步长Simpson积分法 1.算法分析 步长缩小和误

12、差控制方法均与变步长积分雷同，即步长折半、 。基本思想如下： N=1： 求积点 x0=a,x1,x2=b H1=(b-a)/2 S1=f(a)+4*f(x1)+f(b)*H1/3 N=2： 求积点 x0=a,x1,x2,x3,x4=b H2=H1/2 S2=f(a)+4f(x2)+2f(x1)+4f(x3)+f(b)*H2/3 =f(a) +f(b)+2*f(x2) +4*f(x1)+f(x3)*H2/3,5.3 Simpson积分法,N=2n: 关于变步长Simpson积分法的递推计算公式以及程序框图和通用FORTRAN程序将在后面几种积分比较部分提示。这里给出一个通用的FORTRAN程序清

13、单：,(5-8),5.3 Simpson积分法,10 SUBRPUTINE Simpb(a,b,E) 20 h=(b-a)/2 30 s2=0;n=1 40 s0=f(a)+f(b) 45 s1=f(a+h) 50 s=(s0+4.0*s1)*h/3.0 60 n=2*n;h=h/2.0 70 s2=s2+s1 80 s1=0.0 90 x=a+h 100 DO i=1,n 110 s1=s1+f(x),120 x=x+h+h 130 END DO 140 s2n=(s0+2.0*s2+4.0*s1) 150在Simpson积分法中，是用二次三项式（抛物线）去逼近曲线f(x)。而在N-C积分法

14、中，用的是m 次多项式去逼近曲线f(x)。 设在被积区间a, b选取m+1个插值点，用m 次多项式pm去逼近曲线f(x),则： 利用拉格郎日插值多项式,(5-9),5.4 Newton-Cotes积分法,其中: Ck (m)为Cotes系数，可以具体求出起其值： m=1: C0(1) =1/2; C1(1)=1/2 m=2: C0(2)=1/6; C1(2)=4/6; C2(2)=1/6 见下表 注：在实际应用中，一般m9，否则会产生计算不稳定。 上面所谈到的问题，是在整个区间a, b上用Pm去逼近f(x)。下面看看用N-C公式即(5-1)式用在小区间时情况，将a, b分成n个小区间。有：,5

15、.4 Newton-Cotes积分法,5.4 Newton-Cotes积分法,x0=a, x1, x2, , xn-1, xn=b, h=(b-a)/n 当项次m=1时： x0, x1 x1, x2 xn-1,xn,5.4 Newton-Cotes积分法,求和得到 *复化梯形求积公式 同理，当m=2时，在每一个小区间xk, xk+1上增加一个分点（中点），xk+1/2=（xk+ xk+1）/2= xk+ h/2，此时区间a,b上共有2n+1个等距求积节点，即： x0, x0+1/2, x1 , x1+1/2, x2, , xn+1/2, xn 两节点之间的距离（新的步长）h=h/2。因此， m

16、=2时的求积公式为,(5-10),5.4 Newton-Cotes积分法,*复化Simpson求积公式 求和当项次m=4时: 在xk, xk+1上增加三个分点xk+1/4,xk+2/4, xk+3/4 使xk+1/4 = xk+ h/4， xk+2/4 = xk+ 2h/4， xk+3/4 = xk+ 3h/4 ，此时区间a,b上共有4n+1个等距求积节点，两节点之间的距离（新的步长）h=h/4。同理，求和得,(5-12),(5-11),5.5 Romberg积分法,*复化Newton-Cotes求积公式 一、 Romberg积分算法分析 Romberg积分法是用计算精度较低的二分前后两次积分

17、近似值，通过误差补偿法来获得收敛速度较快的一种积分方法。 1、 梯形积分法的误差估计（见式3-3）(线性组合) T*-T2n(T2n-Tn)/3 其中T*为精确值 T* T2n +(T2n-Tn)/3 4*T2n/3-Tn/3 =Sn (5-13)式说明 T2n的误差约等于（T2n-Tn）/3,将这个误差作为一种补偿加到T2n上，可期望得到更精确的Simpson的积分公式 它是梯形积分法二分前后两次积分近似值进行线性组合。,(5-13),(5-14),5.5 Romberg积分法,2、Simpson积分法的误差估计(线性组合) S*-S2n(S2n-Sn)/15 其中S*为精确值 S*S2n

18、+(S2n-Sn)/15 16*S2n/15-Sn/15 =Cn Simpson积分法二分前后两次积分近似值进行线性组合，得到Newton-Cotes的积分公式（精度更高）。 3、Cotes积分公式进行线性加速 C*-C2n(C2n-Cn)/63 其中C*为精确值。 同理 C* C2n +(C2n-Cn)/63 64*C2n/63-Cn/63 =Rn,(5-15),(5-16),(5-17),(5-18),5.5 Romberg积分法,Cotes积分法二分前后两次积分近似值进行线性组合，得到Romberg的积分公式（精度更高）。注意，此时不能再对Romberg的积分进行线性加速。因为此时插值多

19、项式已经为八次，它具有九次插值多项式的精度。依据插值多项式的理论，当插值多项式的次数超过九次时，将出现不稳定情况。因此，到此为止。 Romberg积分法不仅精度高，而且它从最简单的梯形积分开始，逐步加工成收敛速度快的Cotes积分法。 二、 Romberg积分法通用程序框图与通用程序设计 1.算法过程,5.5 Romberg积分法,图5-6 Romberg法子积分图,end subroutine,h=b-a;t1=(f(a)+f(b)*h/2; K=1,s=0; x=a+h/2,s=s+f(x); x=x+h,xb?,t2=t1/2+s*h/2,s2=t2+(t2-t1)/3,k=1?,k=k

20、+1;h=h/2;t1=t2;s1=s2,c2=s2+(s2-s1)/15,k=2?,r2=c2+(c2-c1)/63,c1=c2,k=3?,r1=r2,abs(r2-r1)E?,Y,N,N,N,N,Y,N,Y,Y,subroutine,图5-7 Romberg法子程序框图,5.5 Romberg积分法,10 subroutine romberg(a,b,E) 20 h=b-a;t1=(f(a)+f(b)*h/2 30 k=1 40 s=0;x=a+h/2 50 if(xb) then 60 s=s+f(x);x=x+h/2+h/2 70 endif 80 t2=t1/2+s*h/2;s2=t

21、2+(t2-t1)/3 90 if(k=1) then 100 k=k+1;h=h/2;t1=t2;s1=s2 110 goto 40 115 endif 120 c2=s2+(s2-s1)/15 130 if(k=2) then 140 c1=c2,150 goto 100 160 endif 170 r2=c2+(c2-c1)/63 180 if(k=3) then 190 r1=r2 200 goto 140 210 endif 220 if(abs(r2-r1)E) then 230 goto 190 240 endif 250 end subroutine romberg,3.程序语

22、言,5.5 Romberg积分法,读懂本程序的关键在于看懂控制变量k(步长折减次数）。当k=1时，程序刚开始，程序只求出t2和s2(实际上s1)；当k=2时，步长折二分，t1=t2, s1=s2,且程序又求出新的t2(实际上为t4)和s2，同时求出c2，并且令c1=c2；当k=3时，步长再折二分，t1=t2, s1=s2, 同时又求出新的t2(实际上为t8)和s2(实际上为s4)。求出新的c2,且执行c1=c2,求出r2,且执行r1=r2, c1=c2；当k=4时，步长再折二分， t1=t2, s1=s2, 同时又求出新的t2(实际上为t16)和s2(实际上为s8)以及c2(实际上为c4)，新

23、的r2。此后，只要当k3，程序每次都执行判断abs(r2-r1)是否满足精度要求，若满足，则结束，否则执行r2=r1和c1=c2，且重复k=4（只是k值不同）的过程，直至满足精度要求。,5.6 Gauss积分法,前面介绍的几种积分法都是在等分区间情况下进行讨论的，而Gauss积分法则是适用于不等分区间的一种非常有效的数值积分方法。 一、高斯点及高斯公式简介 首先介绍一下代数精度的概念。一个求积公式，如果对于任何一个不超过m次多项式fm(x)均能使其余项积分为零，即Rf=0，但对于m+1次多项式不能能使其余项积分为零，即Rf0,则称此求积公式具有m次代数精度。 内插求积公式可以表示为,(5-19

24、),5.6 Gauss积分法,在等距节点的情况下，我们可以验证它具有n次代数精度。假定若节点数目一定，我们能否通过适当地选择节点的位置，使积分公式具有更高的代数精度呢？Gauss首先指出，若我们选择的接点x0,x1,x2, xn所构成的多项式 与任意次数不超过n次的多项式p(x)都正交，即 则因 为次数不超过2n+1次的多项式，故求积公式(5-19)式具有2n+1次代数精度。 考虑积分,5.6 Gauss积分法,的上下限a,b是变动的，通过积分变换可将上下限变为固定的值 1和1。事实上，令 ，可以将积分区间a,b变换为-1,1，即 因此，不失一般性，可考察区间-1,1上的Gauss积分公式,(

25、5-20),5.6 Gauss积分法,常用的Gauss积分公式如下： （1）Gauss-Legendre积分公式 （4）Gauss-Hermite积分公式 （2）Gauss-Chebyshev积分公式 （3）Gauss-Laguerre积分公式,5.6 Gauss积分法,表Gauss-Legendre求积公式的xi与Ai,5.6 Gauss积分法,表2 Gauss-Laguerre求积公式的xi与Ai,5.6 Gauss积分法,表3 Gauss-Hermite求积公式的xi与Ai,5.6 Gauss积分法,例用阶数n=2的Gauss求积公式计算以下积分的近似值 解：由于n=2，在表中查出Gau

26、ss求积点坐标和权系数,5.6 Gauss积分法,subroutine,xm=(b+a)/2; xr= (b-a)/2; s=0,DO i=1,n,dx=xr*x(i),Comp f(xm+dx), f(xm-dx),s=s+w(i)*(f(xm+dx)+f(xm-dx),END DO i,s=xr*s,END subroutine,INPUT x(i), w(i),图5-8 Gauss积分子程序框图,5.6 Gauss积分法,10 subroutine gausint(a,b,s) 20 dimension x(i),w(i) 30 data w/0.2955242,0.2692667 40

27、 xr=(b-a)/2 100 s=0 110 do i=1,5 120 dx=xr*x(i),130 s=s+w(i)*(f(xm+dx)+ 140 x2=x+h,DO x=a,b,h,comp f(x1), f(x2),f(x),END DO x,END subroutine,ddf1=(f(x2)-f(x)/h,ddf2=(f(x)-f(x1)/h,ddf3=(f(x2)-f(x1)/2/h,图5-9 差商数值微分子程序框图,5.7 数值微分法,10 subroutine diffder(a,b,h) 20 do x=a,b,h 30 x1=x-h;x2=x+h 40 ddf1=(f(x

28、2)-f(x)/h 50 ddf2=(f(x)-f(x1)/h 60 ddf3=(f(x2)-f(x1)/h 70 end do 80 End subroutine diffder,5.7 数值微分法,二、插值法求数值微分 1、算法分析 用插值公式求数值微分有两点公式、三点公式、四点公式和五点公式等等，下面介绍常用的三点公式。 已知f(x)在等距节点 的函数值 ，且函数f(x)的二阶导数存在，作三点二次插值函数F(x)并对F(x)求导得三点数值导数公式如下：,(5-26),(5-27),5.7 数值微分法,即用这种方法求数值微分时，节点数必须是3的倍数。 要特别强调的是，当插值多项式pn(x)

29、收敛于f(x)时， 不一定能保正也收敛于 ，而且当节点间距h缩小时，虽然截断误差逐步减小，但舍入误差却可能会增大，因此减小步长h不一定能提高计算精度。 2、程序框图与程序设计,(5-28),5.7 数值微分法,subroutine,x0=x;x1=x+h; x2=x+2h,do x=a,b,3h,comp f(x0), f(x1),f(x2),END DO x,END subroutine,idf1=(4f(x1)-3f(x0)-f(x2)/2h,idf2=(f(x2)-f(x0)/2h,图5-10 插值数值微分子程序框图,idf3=( f(x0)-4f(x1)+3f(x2)/2h,5.7 数值微分法,10 subroutine intpder(a,b,h) 20 do x=a,b,3*h 30