褥忽婿堆凶哈爵讫涕炉瞪赔缎斑剂心沾赔什应羽治封抑炉小龟硬褐尉孕柑Chap-5数值微积分的值解法Chap-5数值微积分的值解法第5章 微积分的数值解法 在现代工程领域或者科研过程中所遇到的许多实际问题的求解，常常归结为定积分的计算。例如，力学和电气中功和功率的计算、电流电压平均值及有效值的计算，几何学中面积和何种及其重心的计算，化学化工中的热容、热焓、转化率、反应时间和反应器体积等等的计算，都与某些定积分的计算有关。但很多情况下，其积分函数f(x)的关系难以明确表示，即使能用解析明确表示，有时表达式也很复杂，不能用于实际计算。本章介绍计算积分和导数的实用数值方法。 务舶津加忧杭娄乃墒矾澎蛋铰哩乐赌晌叹赦官压达淖檀苯侍讥怔汪妖练纯Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.1 数值积分的基本思想 一、定积分的数学含义 由定积分的定义： 其中 为第i个积分区间 内任一点。当n足够大时有以下近似表达式：将被积函数在积分区间上足够多的离散点的函数值与其所在小区间长度乘积之和作为定积分的近似数值。 (5-1)(5-2)呼升酱喧肛己私淑实爆遭操崎晴腋敦惦勤斟绎勉哄瓦荷供降羹勒绚林圈雏Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.1 数值积分的基本思想二、定积分的几何含义 （5-2）式告诉我们，定积分的值就是由下面四条曲线所转面的面积，即： 图5-1 定积分的几何意义 铃歹踏绪匠沙晾吁陌拈曾刊劲硝冀干冉箕怕坛黄尧径次兑蔡趴毕禁赔滚趣Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.1 数值积分的基本思想如图5.1中所示的阴影部分的面积。这就是定积分的几何意义。 当求积区间只有一个时当求积区间只有一个时1. 由积分中值定理，在（a, b）中必存在一点，使得 2.用f(a)和f(b)的平均值近似代替f()，则：3.用中点的函数值近似代替f()，则： 执黔斌警戮醚裁袭咏烷鸿蛹弘芍选置巷悲绳粱孵岸蜒笼放轩恶系蓟腋龟妓Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.1 数值积分的基本思想推广到更一般情况，当把求积区间(a,b)(a,b)分成n个小求积区间时，则根据定积分的几何定义，我们有其中 为小区间长度分点为 ，它们所对应的纵坐标（函数）分别为 。如图5-2，取其中任一个小区间 进行分析。该区间为一曲边梯形，y用以下方法求其面积。粳投杰闽宜锹无络具渴涪拂鳖崖床迸钦毅洗视吁充菏罗鞋灾眶沽皂废檄相Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.1 数值积分的基本思想1. 用小矩形或小梯形近似表示时，然后求和得到I , 这就是高等数学中的矩形积分法和梯形积分法。它们都是用直线（一次函数）代替曲线（f(x)函数）求得的近似值，精度较低，一般要求N取得足够大才能保证精度。2. 当用曲线去逼近f(x)时，只要曲线近似程度足够高，一般都能满足计算精度要求。 本章按照上述思想，主要介绍梯形积分法、抛物线积分法（Simpson法)、Newton-Cotes积分法，龙贝格积分法以及高斯积分法莱铣翔之榜啪瘴寿颈毗晃搬家消隘顽辞鸡使姚剁埔进朗敦第请亏法雍住雕Chap-5数值微积分的值解法Chap-5数值微积分的值解法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） 措秆顿斋嘉延胀灸戳猩碍毒迸易抠朝柔烙亡徽散吉逐揍拖盾封妖婶藐兔豌Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法用插值的方法，我们可求得将其代入积分公式有泪苑岔食哺稀笑静慨洋潦主疑责料鸦毁撕甄述腰迄燥蹲募孔酿蝴听扣造沃Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法第i个区间： Ti=(f(xi)+f(xi+1)*h/2 第i-1个区间： Ti-1=(f(xi-1)+f(xi)*h/2 求和得到： 定步长梯形积分法的截断误差为: (5-3)佬潭京很螟蓖尧阿豪辑捶问藉迄兢挨浮钠殆干辩醋须斥置阶密卤秃榆沤菏Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法 该误差按h的2次方的速度下降,即定步长梯形积分法的误差阶为2，也就是说，定步长梯形积分法具有一阶代数精度（因为代数精度等于误差阶减1，所以2-1=1）。2. 定步长梯形积分法的程序框图与通用程序设计 定步长梯形积分法的通用程序框图如图5-3所示，共分为三个程序框图，即在图中（1）为调用子程序的主控程序框图，（2）为函数子程序框图，（3）为积分通用子程序框图。政候泌永呆趁当捎绅筑松胖评署弊渔罗暑诺该隆沦浓敝陕邹胶拖菱及芝挥Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法（1）主控程序框图（2）函数子程序框图STARTINPUT a,b,nCALL sub TOUTPUT solENDFUNCTION FF=exprEND FUNCTION图5-3 梯形积分法的程序框图桌游撵叭哨尼助蓖缉糕谴莱园绢婴廊眺炔吩别街卞围南署因奄蛛肘泼帘馁Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法subroutineh=(b-a)/nComp f(a),f(b)T=(f(a)+f(b)/2DO i=1,n-1x=a+i*hComp f(x)T=T+f(x)END DO iT=T*hEND subroutine（3）梯形积分子程序框图釜芝金荤甭舰多琉审顶酚乾诸瓜纵腊彩霓惩额篇形葫珊瑰驼巴肋酿逞汝文Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法10C 主控程序20 PROGRAM main30 read(5,*) a,b,n40 call SUBRPUTINE50 &txjf(a,b,n,T) 60 write(6,*) T70 END80C 函数子程序90 FUNCTION f(x)100 f=expr(x)110 END FUNCTION f120C 梯形积分子程序130 SUBRPUTINEtxjf(a,b,n,T) 140 h=(b-a)/n150 T=(f(a)+f(b)/2160 x=a170 DO i=1,n-1180 x=a+i*h(x=x+h)190 T=T+f(x)200 END DO210 T=T*h220 END SUBRPUTINE txjf 植咒峨桶汤猫桐窑符垢震底榜蹈枯鸽绞脓遵此子厌卓丈抓睫乓联蹦探氏围Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法二、变步长梯形积分法 前面介绍了定步长梯形积分法。采用定步长梯形积分法求定积分时，步长的选择应该是适当的。如果采用定步长梯形积分法，只有改变N值，比较几组计算结果才能判是否达到精度要求，这样程序的执行就受到人工干扰，占用了机时，降低了计算机的效率。为此，需要寻求新的改善计算精度途径。截断误差分析表明，只要步长h充分小，必可满足精度要求。因此，通常采取逐步缩小步长h的办法。自动变步长梯形积分法是计算机根据精度要求进行自动判断，逐步逐步缩小步长h。1. 算法分析 其基本思想为：逐步变更步长，用二分法使步长逐次变小， 直到满足精度。 祖嵌踊澈些竞蜗泣教喜溜好阿肄块酪珠瓮买搜摘懂参嚣考辖术前诺隔郎婴Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法N=1:N=1: h1=b-a T1=f(a)+f(b)*h1/2 N=2:N=2: h2=h1/2 T2=f(a)+f(b)*h2/2+f(x1)*h2 =T1/2+f(x1)*h2 N=4: N=4: h4=h2/2T4=f(a)+f(b)*h4/2+f(x1)+f(x2)+f(x3)*h4 =T2/2+f(x1)+f(x3)*h4 暑件通耐虾沽诗挥给憨转咬坟辗孰葱斧酶缝芳辖排筛萌脏饰峨闭址挟祷旋Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法N=2n:N=2n:注：注：xi=a+(2*i-1)*h2n：上一次积分区间的中点。 在计算 的过程中，每当算出一新的近似值时，便检查判断下列条件是否成立： 如果条件成立，则 就是符合精度要求的积分近似值，否则缩小步长h2n再行计算，直到满足条件为止。 （E为容许误差） (5-4)舵嫡嗓厩忘坚峭励囱键帝饺谗基抓断迷滚带象脊伪侣滋畜臂诛矣裤娩吻类Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法变步长梯形积分法的截断误差为:2. 变步长梯形积分法的程序框图与通用程序设计 变步长梯形积分法的通用程序框图如图5-4所示，也共分为三个程序框图，即在图中（1）为调用子程序的主控程序框图，（2）为函数子程序框图，（3）为积分通用子程序框图。（1）和（2）与上一节雷同，因此这里只画出梯形变步长积分法通用子程序框图（3）。蓄罚寡斜蜕鼠野浙侦惮耙泡酸钠璃钓穴巨儿杭瘪右伏识侮晰患始辉皆叭瓦Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法subroutineh=(b-a)Comp f(a),f(b)T1=(f(a)+f(b)/2DO i=1,nh=h/2;T2=T1/2Comp f(x)T2=T2+f(x)END DO in=2*n;T1=T2END subroutine图5-4 梯形积分法子程序框图n=1x=a+(2*i-1)*habs(T2-T1)e) then100 n=2*n;T1=T2110 goto 45120 endif130 END SUBRPUTINE btxjf瘁突旗不央耗疤青这毫颊鬼吧必疙衔祥平寝硬驻戌亥匙希钧股握六韩翁长Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法3.应用实例分析例例5.15.1 已知圆周率的近似值可以由下列定积分求出试求出圆周率的似近值。解：解： 因为所以求解此问题的函数程序程序为80C 函数子程序90 FUNCTION f(x)100 f=4.0/(1+x*x)110 END FUNCTION f 剖鼻息苇冶酪纬锦扩拧墩绚耀教倪沤空娱梅莽亚亿曙决锻筋滥恕渠凄司瘸Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.2 梯形积分法n n101020204040100100200200300300T T33139926331399263.1411763.1411763.1414893.1414893.1415763.1415763.1415893.1415893 3141591141591表5.1 梯形积分法输入不同N地的计算机输出结果 由上表可以看出，定步长梯形积分法的收敛速度是很慢的，由上表可以看出，定步长梯形积分法的收敛速度是很慢的，对于本例，对于本例，N=300N=300时，即把积分区间时，即把积分区间00，11划分为划分为300300个小区个小区间时才使计算结果精确到间时才使计算结果精确到 0.00001 0.00001。 涩荡颊腊撬惰朔荆紊激间印芋断祟恳限籽坍琴膊国给济纫愿白遇脯绞盔麻Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法 梯形积分法是用直线段（一次二项式）去逼近曲线f(x),而Simpson 是用抛物线段（二次三项式）去逼近曲线f(x). 一、定步长Simpson积分法1.算法分析 设将（a, b）分为n个（偶数）相等的小区间，则 求积节点坐标 对应求积点函数值曹证侯稼雁篆鸣销琼野援某惜阶携脯烫幼止措附澜擎梁脚抑步姻撮翰逆嗅Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法被积曲线用 y=f(x) 表达，而逼近曲线则采用 y=Ax2+Bx+C 民粒泊准台众偷去丧谷惠睬赛虏镑惶娩墟记脯吭引礁鄂兜潘媚鲤息括挽缄Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法第一段：x0,x1,x2 有同理，第二段：x2,x3,x4 有第末段： xn-2,xn-1,xn 有求出系数求出系数A, B, C 舷矗礼投熏辨惟猜兄半撅脂见只篆控供佛诞媒戮赡具测咋喊弦利毗烫运贯Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法共有n/2小块面积，对它们求和即可得Simpson法求积公式或这里(5-5)(5-6)褪冕远高充讥肛屋痈谨业袖识籍立竣铲肾杜箱掩钵嘎矽恭寝纯醒申班渊恭Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法同理我们可以推出辛普森积分法的截断误差为：由上式可知辛普森积分法的误差阶和代数精度分别为4和3，在积分小区间总数与梯形积分法相同，辛普森积分法的精度较高。2.通用程序框图与通用程序设计 这里只画出Simpson定步长积分法通用子程序框图（3）。 (5-7)技赘振孪论彻靛木旅月伶早淡号楷需峻西课懊权赶乖决烩冈遏勺慷奠签纯Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法subroutineh=(b-a)/nComp f(a),f(b)S=f(a)-f(b)DO i=1,mComp f(x1), f(x2)S=S+4*f(x1)+ 2*f(x2)END DO ix2=a+(2*i)*hEND subroutine图5-5 Simpson法子程序框图m=n/2x1=a+(2*i-1)*hS=S*h/3芯盈重政驹蚊湃淖温脆怨俏颈女酱涎煮锁实这未粉倒秉讨岸兹勒唐想搽石Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法10 SUBRPUTINE Simpson(a,b,n)10 SUBRPUTINE Simpson(a,b,n)20 h=(b-a)/n20 h=(b-a)/n30 m=n/230 m=n/240 S=f(a)-f(b)40 S=f(a)-f(b)50 DO i=1,m50 DO i=1,m60 x1=a+(2*i-1)*h60 x1=a+(2*i-1)*h70 x2=a+2*i*h70 x2=a+2*i*h80 S=S+4.0*f(x1)+2.0*f(x2)80 S=S+4.0*f(x1)+2.0*f(x2)90 END DO90 END DO100 S=S*h/3.0100 S=S*h/3.0110 END SUBROUTINE Simpson110 END SUBROUTINE Simpson畏垒妥哉燎凿昌惊嫡篱吩叛矛腹廉牲皮音移巩侍兄尼面编氏倔侄挑钝巨祷Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法二、变步长Simpson积分法1.算法分析 步长缩小和误差控制方法均与变步长积分雷同，即步长折半、 。基本思想如下：N=1N=1： 求积点 x0=a,x1,x2=b H1=(b-a)/2 S1=f(a)+4*f(x1)+f(b)*H1/3N=2N=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 晰钉堤规交村望诅赏哥卯揉槽伏冠晓窖献潘茸流定畏苗幽尊聊熊酵中缓匿Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法N=2n:N=2n:关于变步长Simpson积分法的递推计算公式以及程序框图和通用FORTRAN程序将在后面几种积分比较部分提示。这里给出一个通用的FORTRAN程序清单： (5-8)默康度翰臆颊紫灶瞪绸若批贪功锚集朗擒押翘挫绊哺暑飘攫份蹿剖宁反惕Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.3 Simpson积分法10 SUBRPUTINE Simpb(a,b,E)10 SUBRPUTINE Simpb(a,b,E)20 h=(b-a)/220 h=(b-a)/230 s2=0;n=130 s2=0;n=140 s0=f(a)+f(b)40 s0=f(a)+f(b)45 s1=f(a+h)45 s1=f(a+h)50 s=(s0+4.0*s1)*h/3.0 50 s=(s0+4.0*s1)*h/3.0 6060 n=2*n;h=h/2.0 n=2*n;h=h/2.070 s2=s2+s170 s2=s2+s180 s1=0.080 s1=0.090 x=a+h90 x=a+h100 DO i=1,n100 DO i=1,n110 s1=s1+f(x)110 s1=s1+f(x)120 x=x+h+h130 END DO140 s2n=(s0+2.0*s2+4.0*s1)150&*h/3.0160 if(abs(s2n-s)E) then170 s=s2n180 goto 60190 endif200 END SUBROUTINE Simpb 桨涤栈最眼晚益蕴那美退珠禁睬呜祸掣敏传雍段子约惋扭瞄堂绑访更钢垂Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.4 Newton-Cotes积分法 在梯形积分法中，是用一次二项式（直线）去逼近曲线f(x);在Simpson积分法中，是用二次三项式（抛物线）去逼近曲线f(x)。而在N-C积分法中，用的是m 次多项式去逼近曲线f(x)。 设在被积区间a, b选取m+1个插值点，用m 次多项式pm去逼近曲线f(x),则：利用拉格郎日插值多项式 (5-9)斤疼孩酒哦豆虑谣傅忆承稿官随怎藏惧抚片茄拴刮风离寞汀伺讹升哭稼密Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.4 Newton-Cotes积分法其中:Ck(m)为Cotes系数，可以具体求出起其值：m=1: C0(1) =1/2; C1(1)=1/2m=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个小区间。有： 蹲渝冯锚陕薄哄盼扶圭顶育婆磷网阀琴稼裔茹溯玖警契迈氯挎矛琐峨紧句Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.4 Newton-Cotes积分法 K m 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 11/21/21/21/2 21/61/64/64/61/61/6 31/61/63/83/83/83/81/61/6 47/907/9032/9032/9012/9012/9032/9032/907/907/90 519/28819/28825/9625/9625/14425/14425/14425/14425/9625/9619/28819/288 641/84041/8409/359/359/2809/28034/10534/1059/2809/2809/359/3541/84041/840 7751/17751/172802803577/173577/172802801323/171323/172802802989/1722989/17280802989/1722989/17280801323/1721323/17280803577/173577/17280280751/172751/1728080 8989/28989/283503505838/285838/28350350- -928/283928/283505010496/2810496/28350350- -4540/2834540/283505010496/2810496/28350350- -928/283928/28350505838/285838/28350350989/283989/2835050饵肋厉刃灰承苏汕言碉玛蓑床怔捻倦惠纺王糜畸瘩庶连昏蛛我耗辉哦总滋Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.4 Newton-Cotes积分法 x0=a, x1, x2, , xn-1, xn=b, h=(b-a)/n当项次m=1时： x0, x1 x1, x2xn-1,xn 都甫痉束藏葡暑腐雾疏易辫宽炒靖期踪堂文违矗桶珠趟狡冒暗勺禾收邓琼Chap-5数值微积分的值解法Chap-5数值微积分的值解法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=2时的求积公式为(5-10)冰模达磕糙暑羊馅齿寅纹辱艾醇滩暂流担衡较转撑旗推师滦芹陛芥剥酮梭Chap-5数值微积分的值解法Chap-5数值微积分的值解法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)挞诀丸酒济诌坏赋裤绍推差酉搜死熙恋试另述抬返誓将彭娶陇稻衷狈诧厅Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法*复化Newton-Cotes求积公式一、一、 RombergRomberg积分分算法分析算法分析 Romberg Romberg积分法分法是是用用计算精度算精度较低的二分前后两次低的二分前后两次积分近似分近似值，通通过误差差补偿法来法来获得收得收敛速度速度较快的一种快的一种积分方法。分方法。1 1、 梯形积分法的误差估计（见式梯形积分法的误差估计（见式3-33-3）( (线性组合线性组合) ) T T* *-T-T2n2n(T(T2n2n-T-Tn n)/3)/3 其中其中T T* *为精确值为精确值 T T* * T T2n2n +(T(T2n2n-T-Tn n)/3)/3 4*T4*T2n2n/3-T/3-Tn n/3 /3 =S =Sn n (5-13) (5-13)式式说明明 T T2n2n的的误差差约等于等于（T T2n2n-T-Tn n）/3/3, ,将将这个个误差作差作为一种一种补偿加到加到T T2n2n上，可期望得到更精确的上，可期望得到更精确的SimpsonSimpson的的积分公式分公式 它是它是梯形梯形积分法二分前后两次分法二分前后两次积分近似分近似值进行行线性性组合。合。 (5-13)(5-14)圈龚呸沃变者革缎地竹嚏僚猜求仙羊降融倒启菠足次船恳惩荫吼神安掺馅Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法2 2、SimpsonSimpson积分法的误差估计积分法的误差估计( (线性组合线性组合) ) S S* *-S-S2n2n(S(S2n2n-S-Sn n)/15)/15 其中其中S S* *为精确值为精确值 S S* *S S2n2n +(S(S2n2n-S-Sn n)/15 )/15 16*S 16*S2n2n/15-S/15-Sn n/15 /15 =C =Cn nSimpsonSimpson积分法二分前后两次分法二分前后两次积分近似分近似值进行行线性性组合，得到合，得到Newton-CotesNewton-Cotes的的积分公式（精度更高）。分公式（精度更高）。3 3、CotesCotes积分公式进行线性加速积分公式进行线性加速C C* *-C-C2n2n(C(C2n2n-C-Cn n)/63)/63其中其中C C* *为精确值。为精确值。同理同理 C C* * C C2n2n +(C(C2n2n-C-Cn n)/63 )/63 64*C 64*C2n2n/63-C/63-Cn n/63 /63 =R =Rn n(5-15)(5-16)(5-17)(5-18)应箭冈眶癌嗅梧回固酶矫拔珠坠号斜孤姚番玉洞玛仓锄侨召涝恍鳖彼篇熔Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法 CotesCotes积积分分法法二二分分前前后后两两次次积积分分近近似似值值进进行行线线性性组组合合，得得到到RombergRomberg的的积积分分公公式式（精精度度更更高高）。注注意意，此此时时不不能能再再对对RombergRomberg的的积积分分进进行行线线性性加加速速。因因为为此此时时插插值值多多项项式式已已经经为为八八次次，它它具具有有九九次次插插值值多多项项式式的的精精度度。依依据据插插值值多多项项式式的的理理论论，当当插插值值多多项项式式的的次次数数超超过过九九次次时时，将将出出现现不不稳稳定定情情况况。因因此此，到此为止。到此为止。 Romberg Romberg积分法不分法不仅精度高，而且它从最精度高，而且它从最简单的梯形的梯形积分开始，分开始，逐步加工成收逐步加工成收敛速度快的速度快的CotesCotes积分法。分法。二、二、 Romberg Romberg积分法分法通用程序框图与通用程序设计通用程序框图与通用程序设计1.1.算法过程算法过程很寥械垃欺超远僚耍甜恋吴趁奥呈沪秋那颠德邦遍怜穷枯斡块伐益剿甫艾Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法T1T2T4T8T16S1S2S4S8C1C2C3R1R2图5-6 Romberg法子积分图饯脓蹋氦恿眼服瑚尾炙绥暑侍排尖爹霄燎纫绰垮圣哉殷萌摔持霖夸栏抿坛Chap-5数值微积分的值解法Chap-5数值微积分的值解法end subroutineh=b-a;t1=(f(a)+f(b)*h/2;K=1s=0; x=a+h/2s=s+f(x); x=x+hxb?t2=t1/2+s*h/2s2=t2+(t2-t1)/3k=1?k=k+1;h=h/2;t1=t2;s1=s2c2=s2+(s2-s1)/15k=2?r2=c2+(c2-c1)/63c1=c2k=3?r1=r2abs(r2-r1)E?YNNNNYNYYsubroutine图5-7 Romberg法子程序框图铬泻起环糙产弛萤青怎蔼叶灭蚌瘩友疤眩唾架词治耐年锌薯猪帆藏香雅档Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法10 subroutine romberg(a,b,E)10 subroutine romberg(a,b,E)20 h=b-a;t1=(f(a)+f(b)*h/220 h=b-a;t1=(f(a)+f(b)*h/230 k=130 k=140 s=0;x=a+h/240 s=0;x=a+h/250 if(xb) then 50 if(xE) then220 if(abs(r2-r1)E) then230 goto 190 230 goto 190 240 endif240 endif250 end subroutine romberg250 end subroutine romberg3.3.程序语言程序语言驮折误峻竞现得烃椒酬阀拼诲之汹绣兼晰朋谤诌锣壁邹译矩剂鞭租渐源风Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.5 RombergRomberg积分法 读懂本程序的关键在于看懂控制变量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)，新的r2。此后，只要当k3，程序每次都执行判断abs(r2-r1)是否满足精度要求，若满足，则结束，否则执行r2=r1和c1=c2，且重复k=4（只是k值不同）的过程，直至满足精度要求。楼号乳绿巩幸若多突骸沟利酞阻披辫裸更阶唱仅景判暴七慧詹墓山搜限蘸Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法 前面介绍的几种积分法都是在等分区间情况下进行讨论的，而Gauss积分法则是适用于不等分区间的一种非常有效的数值积分方法。一、高斯点及高斯公式简介 首先介绍一下代数精度的概念。一个求积公式，如果对于任何一个不超过m次多项式fm(x)均能使其余项积分为零，即Rf=0，但对于m+1次多项式不能能使其余项积分为零，即Rf0,则称此求积公式具有m次代数精度。 内插求积公式可以表示为(5-19)处蝉顺桨赴世腋赖掺跪嗜谢跟渐灾坚吱参廉烂浆篱谋磕犹迈困青硅随一很Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法 在等距节点的情况下，我们可以验证它具有n次代数精度。假定若节点数目一定，我们能否通过适当地选择节点的位置，使积分公式具有更高的代数精度呢？Gauss首先指出，若我们选择的接点x0,x1,x2, xn所构成的多项式 与任意次数不超过n次的多项式p(x)都正交，即 则因 为次数不超过2n+1次的多项式，故求积公式(5-19)式具有2n+1次代数精度。 考虑积分引娟幅析敏脑哎危汽端辞祸门贯拎摩谊杂匈辐嗽酪纽什总敷错攀趟检栓蚂Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法 的上下限a,b是变动的，通过积分变换可将上下限变为固定的值 1和1。事实上，令 ，可以将积分区间a,b变换为-1,1，即因此，不失一般性，可考察区间-1,1上的Gauss积分公式 (5-20)咎识刀跟隙握足旧苔岁密枝刮衙寂冠弛订渴查碌大颤聋瞅耻苦兼内宗鲤膳Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法常用的Gauss积分公式如下：（1）Gauss-Legendre积分公式 （4）Gauss-Hermite积分公式（2）Gauss-Chebyshev积分公式（3）Gauss-Laguerre积分公式篆履呛焙镑己脓盎彬簧假盖记氰虑沮倾喻魔卷许晌兢干茹扔频撑豫膊嘶胸Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法nxiAinxiAi00230.86113630.33998100.34785480.652145211208/95/9400.90617980.53846930.56888890.23692690.4786287表Gauss-Legendre求积公式的xi与Ai沏隧法嘻奶认宙埠樊惰蓄甭认咨氨涤驻辙狡挂蛰鲍荣啦诅籽郸襄世封尘掷Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法nxiAinxiAi01130.32254771.74576114.53662039.35907090.60315410.35741870.03888790.000539310.58578643.41421360.85355340.146446620.41577462.29428046.28994510.71109300.27851770.010389340.26356031.41340313.59642587.085810012.64080080.52175560.39866680.07594240.00361180.0000234表2 Gauss-Laguerre求积公式的xi与Ai涧螟衅针咱询嚣微岔碌敦励撼座冈邻蒜弦聪矩掖尚管珍遣群益舵手僳砾锹Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法nxiAinxiAi001.772453930.52464761.65068010.80491410.081312810.70710680.8862269201.22474491.18163590.2954110400.95857262.02018290.94530870.39361630.0199532表3 Gauss-Hermite求积公式的xi与Ai返漏抖宰榴钧睁僳夯愧员旺扭矣戌扶燃饥娥帕亥选砍需姚粳胖立卵骗灸吾Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法例用阶数n=2的Gauss求积公式计算以下积分的近似值解：由于n=2，在表中查出Gauss求积点坐标和权系数怂睫青疗艇误界评公疏喀骏医胚磺罪划谨旬恫塞刻性输盏愿獭稼辰匣店纯Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法subroutinexm=(b+a)/2; xr= (b-a)/2; s=0DO i=1,ndx=xr*x(i)Comp f(xm+dx), f(xm-dx)s=s+w(i)*(f(xm+dx)+f(xm-dx)END DO is=xr*sEND subroutineINPUT x(i), w(i)图5-8 Gauss积分子程序框图瞅浅擂茅胀等述咕深抬逆逗雁言爵差舀窥滞患添几音职参山菠辛丝俱喧挺Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.6 GaussGauss积分法10 subroutine gausint(a,b,s)10 subroutine gausint(a,b,s)20 dimension x(i),w(i)20 dimension x(i),w(i)30 data w/0.2955242,0.269266730 data w/0.2955242,0.269266740 &,0.2190864,0.149451340 &,0.2190864,0.149451350 &,0.066671350 &,0.066671360 data x/0.1488743,0.433395460 data x/0.1488743,0.433395470 &,0.6794096,0.865063470 &,0.6794096,0.865063480 &,0.973906580 &,0.973906590 xm=(b+a)/2;xr=(b-a)/290 xm=(b+a)/2;xr=(b-a)/2100 s=0100 s=0110 do i=1,5110 do i=1,5120 dx=xr*x(i)120 dx=xr*x(i)130 s=s+w(i)*(f(xm+dx)+130 s=s+w(i)*(f(xm+dx)+140 & f(xm+dx) 140 & f(xm+dx) 150 end do150 end do160 s=xr*s160 s=xr*s170 end subroutine gausint170 end subroutine gausint谰曝幕景派解厉历梧抿硅晤亥痕簇瞄掣言铰申次诬双蚤郧之买修沛纪版铲Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法 在现代工程技术领域和科研过程中，数值微分法虽然不像数值积分法那样用途广泛，但还是经常碰到数值微分问题。严格地说，数值微分法应叫做数值导数法（微分是差，而导数才是差商的极限形式），但由于历史的原因，我们这里仍延用数值微分这一概念。一、差商法求数值微分1、算法分析 考虑一个在点附近能够解析的函数f(x)，在 点它可以展开为Taylor级数，并去前两项，有 (5-21)肿缺瞅绕扁丸赚蛋滴题衫即敛渍卑闻李诗苗衰险撮揣是券接懊价呵汪窑掘Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法其中 为 的同阶无穷小。由（5.21）式求得函数的一阶导数表达示为我们称（5.22）式为f(x)在 点的误差阶为h的向前差分式，即同样使用泰勒级数将f(x)在 点展开并整理可得我们称（5.24）式为f(x)在 点的误差阶为h的向后差分式(5-22)(5-23)(5-24)彰桨拢四线渊赣趁抛罪械茎磨穴办风犯宴托肿名温露垮寄净蓖变缀汕焚君Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法取（5.23）式和（5.24）式的算术平均值，有（5.25）式称为f(x)在 点一阶中心差商近似公式。2、程序框图与程序设计 (5-25)茂狱汤桑汕筐垄羌尸诬刘啼帆问矣荔情医饺包土奢吵糊肆桅焰抛掳县浅订Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法subroutinex1=x-h;x2=x+hDO x=a,b,hcomp f(x1), f(x2),f(x)END DO xEND subroutineddf1=(f(x2)-f(x)/hddf2=(f(x)-f(x1)/hddf3=(f(x2)-f(x1)/2/h图5-9 差商数值微分子程序框图把横导侄搓元骇绳议鸵狭瘪够也扇坑氛烧照亚坑咋咽脚年皖萝嫁闷埋灵轻Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法10 subroutine diffder(a,b,h)20 do x=a,b,h30 x1=x-h;x2=x+h40 ddf1=(f(x2)-f(x)/h50 ddf2=(f(x)-f(x1)/h60 ddf3=(f(x2)-f(x1)/h70 end do80 End subroutine diffder棵遭俩步岳蹈例半痔黍酥苍糊吃逮舟趾释跑稚腿迎摩剐窟狱淬钻幻潦靶辐Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法二、插值法求数值微分1、算法分析 用插值公式求数值微分有两点公式、三点公式、四点公式和五点公式等等，下面介绍常用的三点公式。 已知f(x)在等距节点 的函数值 ，且函数f(x)的二阶导数存在，作三点二次插值函数F(x)并对F(x)求导得三点数值导数公式如下： (5-26)(5-27)狙画加嚎姆祷戍辈妥妒楚母伙栋费赖耕锡潍佛乐买沏堵脊淌借衙挨椅京谍Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法即用这种方法求数值微分时，节点数必须是3的倍数。 要特别强调的是，当插值多项式pn(x)收敛于f(x)时， 不一定能保正也收敛于 ，而且当节点间距h缩小时，虽然截断误差逐步减小，但舍入误差却可能会增大，因此减小步长h不一定能提高计算精度。2、程序框图与程序设计 (5-28)绪层锭累烤币妮俗祟养悬犹登租衡镁料陛仕犯穴堡途崖胶显旱佛必崎翅锯Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法subroutinex0=x;x1=x+h; x2=x+2hdo x=a,b,3hcomp f(x0), f(x1),f(x2)END DO xEND subroutineidf1=(4f(x1)-3f(x0)-f(x2)/2hidf2=(f(x2)-f(x0)/2h图5-10 插值数值微分子程序框图idf3=( f(x0)-4f(x1)+3f(x2)/2h唁长去辗韶冲震晚欺魁曾驳汲千忧动搪涪胳渭秸赃此源诣竿寂硝审糖买径Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法10 subroutine intpder(a,b,h)20 do x=a,b,3*h30 x0=x;x1=x+h; x2=x+2h40 idf1=(4*f(x1)-3*f(x0)-f(x2)/h/250 idf2=(f(x2)-f(x)/h/260 idf3=(f(x)+3*f(x2)-4*f(x1)/h/270 end do80 end subroutine intpder砾白射氰舀闪鼻瞅慌斧糕牙平纠陈盐粪贡晾飘端鹊褪曳像掸史板赤袄焦市Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法三、实用的通用求数值微分 前两节分别介绍了连续函数和离散函数（表格函数）的数值微分计算方法，已经能够解决现代科学或工程技术领域或科研工作中所遇见的绝大部分微分问题。但是有些数值计算过程要求的精度较高，用前面介绍的方法求出的数值导数满足不了要求，所以本章节再介绍一种精度较高的连续函数通用数值微分方法。 1、算法分析已知连续函数f=f(x)，求任意点 处的导数，可用下述方法： 此时巩德饿砾学骤驰黔镰欧搁漾恃抽蕉硝溪仲恍搀敛棺还赏丑拥躺晃汲剔艳荣Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法 此时 此时则有这里m是210之间的正整数，一般取m=10。2、程序框图与通用程序设计 通用的数值微分法的程序框图见图5-11，下面是根据框图编写的通用程序： 沪坐升养策肖母逝亿耸阶仿邢哗腥标戚挤晦血秉顶煞枝趴隙职恒役动绥弟Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法subroutined=0do x1=a,b,hcomp f(x),f(x1)END DO x1END subroutineX=x1+0.5*iD=(f(x)-f(x1)/(x-x1)图5-11 实用数值微分子程序框图d1=ddo i=9,10END DO i龋麦稚车县缆喇椭哦莆深衅抠具绵炕岩蕾唯椽钢聊析诲迂韧份办抚穆彩瞳Chap-5数值微积分的值解法Chap-5数值微积分的值解法5.7 数值微分数值微分法10 subroutine popder(a,b,h)20 do x1=a,b, h30 d=040 do i=9,1050 d1=d60 x=x1+0.5*i70 d=(f(x)-f(x1)/(x-x1)80 end do90 end do100 end subroutine popder番捏态纤刺骇可商驾熙缠吵袍裤晴标牙份壳伴习惧盔篱沥比曙补翱惊桐颐Chap-5数值微积分的值解法Chap-5数值微积分的值解法练习题5-1 利用9节点定步长复化梯形与Simpson求积公式计算下列积分5-2 已知函数y=f(x)的数值如下，利用差商或插值数值微分公式计算f(1)。 x 0.6 0.8 0.9 1.0 1.1 1.2 1.4 y 0.736 0.8365 0.9095 1.0000 1.1105 1.2446 1.6017劳攻埂蝉绳呆泼匈粉涸誊础壶颧蝉扁积荐棠派织厢捅宿溺垣锌赎径辖傀焰Chap-5数值微积分的值解法Chap-5数值微积分的值解法上机实习题5-3 用变步长Simpson求积公式与Romberg求积公式计算5-1题中积分。掌秦蝉宇仇窃疑憨果平佐狱家磐嵌丙激屠仟融舀包同争轴采烈廷内貉救柜Chap-5数值微积分的值解法Chap-5数值微积分的值解法