第二章数值微分和数值积分

1、数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS第二章 数值微分和数值积分数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS数值微分hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0001. 函数f(x)以离散点列给出时，而要求我们给出导数值，2. 函数f(x)过于复杂这两种情况都要求我们用数值的方法求函数的导数值微积分中，关于导数的定义如下：自

2、然，而又简单的方法就是，取极限的近似值，即差商数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICShxfhxfxf)()()( 000向前差商x0 x0+h数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS由Taylor展开hxxfhxhfxfhxf002000),( ! 2)( )()(因此，有误差)()( ! 2)()()( )(000hOfhhxfhxfxfxR数 学 系University

3、of Science and Technology of ChinaDEPARTMENT OF MATHEMATICShhxfxfxf)()()( 000向后差商x0-hx0数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS由Taylor展开hxxfhxhfxfhxf002000),( ! 2)( )()(因此，有误差)()( ! 2)()()( )(000hOfhhhxfxfxfxR数 学 系University of Science and Technology of ChinaDEPART

4、MENT OF MATHEMATICShhxfhxfxf2)()()( 000中心差商x0-hx0 x0+h数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS由Taylor展开23000010102300002020()()()()( ),2!3!()()()()(),2!3!hhf xhf xhfxfxfxxhhhf xhf xhfxfxfxhx因此，有误差)()( 6)( )( 12 2)()()( )(22212000hOfhffhhhxfhxfxfxR数 学 系University of

5、 Science and Technology of ChinaDEPARTMENT OF MATHEMATICSf(x)=exp(x)hf(1.15)R(x)hf(1.15)R(x)0.103.1630-0.00480.053.1590-0.00080.093.1622-0.00400.043.1588-0.00060.083.1613-0.00310.033.1583-0.00010.073.1607-0.00250.023.1575-0.00070.063.1600-0.00180.013.1550-0.0032例：数 学 系University of Science and Techn

6、ology of ChinaDEPARTMENT OF MATHEMATICS由误差表达式，h越小，误差越小，但同时舍入误差增大，所以，有个最佳步长我们可以用事后误差估计的方法来确定设D(h),D(h/2) 分别为步长为h,h/2 的差商公式。则)2()(hDhD时的步长h/2 就是合适的步长( )( )( )( )( /2)( /2)fxD hO hfxD hO h( )( )( )2( )( /2)( /2)fxD hO hfxD hO h( )( )2( )2 ( /2)fxD hfxD h( )( /2)( )( /2)fxD hD hD h数 学 系University of Sci

7、ence and Technology of ChinaDEPARTMENT OF MATHEMATICS 插值是建立逼近函数的手段，用以研究原函数的性质。因此，可以用插值函数的导数近似为原函数的导数)()()()(xLxfknk误差插值型数值微分用Taylor展开分析数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS给定点列20)(,(iiixfx且hxxxx0112，求)( ),( ),( 012xfxfxf解：)(2)()()()(2)()(2210122002212xfhxxxxxfhx

8、xxxxfhxxxxxL)(2)()()()(2)()(2210122002212xfhxxxxxfhxxxxxfhxxxxxL例：数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS20200121() ()3 ()4 ()()( )23hfxLxf xf xf xfh2121021() ()()()( )26hfxLxf xf xfh22220121() ()()4 ()3 ()( )23hfxLxf xf xf xfhTaylor展开分析，可以知道，它们都是)(2hO称为三点公式三点公式02

9、001221() ()()2 ()()fxLxf xf xf xh11201221() ()()2 ()()fxLxf xf xf xh22201221() ()()2 ()()fxLxf xf xf xh误差？数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS数值积分)()()(aFbFdxxfba关于积分，有Newton-Leibniz 公式但是，在很多情况下，还是要数值积分：1、函数有离散数据组成2、F(x)求不出3、F(x)非常复杂定义数值积分如下：是离散点上的函数值的线性组合)()(0

10、iniinxfafI称为积分系数积分系数，与f(x)无关，与积分区间和积分点有关数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)()()(0afabfI例：)2()()()(2)()(00bafabfIafabfI数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)()(0iniinxfafI为数值积分，badxxffI)()(为积分，则称数值积分有k阶代数精度阶代数精度是指：)()(;, 0

11、),()(11kkniinxIxIkixIxI问题：如果判断好坏？代数精度代数精度 对任意次数不高于k次的多项式f(x)，数值积分没有误差数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS用插值函数的积分，作为数值积分)()()()()()(00inibaiibaniibannxfdxxldxxfxldxxLfI ia代数精度代数精度由Lagrange插值的误差表达式，)()!1()()()1(xnfxRnnn，有dxxnfdxxRfIfIbannbann)()!1()()()()()1(可以看

12、出，至少n 阶代数精度nkxxfxfkn,)(, 0)()1(插值型数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS121111122101010nabababaaaxxxxxxnnnnnnnnVandermonde行列式使用尽可能高的代数精度mixIxIiin, 0),()( niia0已知求系数所以，如果mn，则系数唯一前面得到的系数是最好的吗？数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATI

13、CS若数值积分至少n阶代数精度，则系数唯一( ),0,biiaal x dx in误差误差dxxnxfdxxRfIfIbannbann)()!1()()()()()1(数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSNewton-Cotes 积分若节点可以自由选取，则，一个自然的办法就是取等距节点。对区间做等距分割。该数值积分称为Newton-Cotes积分数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHE

14、MATICSniihaxnabhi,0,dtntititttnininhdthinintititttdxxlaninninthaxbaii00)() 1)(1() 1()!( !) 1( ) 1()!( !)() 1)(1() 1()()()(niiCaba设节点步长)( niC(b-a)与步长h无关，可以预先求出数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSN1时2121) 1(10)1(110)1(0dttCdttC)(21)()(21)()(1bfabafabfI梯形公式数 学 系Un

15、iversity of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSN2 时61) 1(4164)2(2161)2)(1(4120)2(220)2(120)2(0dtttCdtttCdtttC)(61)()2(64)()(61)()(2bfababfabafabfISimpson公式数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS1、梯形公式1( )( )()()2!bxafE fx a x bdx此处用了积分中值定理误差误差

16、3( )()()()( )2!12bafb ax a x bdxf 数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS2、Simpson 公式)()()()()()()(332fSPSPIfIfSfIfE 注意到，Simpson 公式有3 阶代数精度，因此为了对误差有更精确地估计，我们用3 次多项式估计误差)2()2(),()(),()(333bafbaPbfbPafaP为0)(2880)()(2)(! 4)( )(2)(! 4)( )4(52)4(2)4(fabdxbxbaxaxfdxbxba

17、xaxxfbaba)2( )2( ,3bafbaP数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS2)4(32)(! 4)()(! 3)( )(2)( )( )()(cxxfcxcfcxcfxxcfcfxf数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS一般的有kndxxxKfnKfEbannnnn2, 0)(),()!2()()2(12, 0)(),()!1()()1(kndxxKfnKfE

18、bannnnn因此，N-C积分，对偶数有n+1 阶代数精度，而奇数为n 阶代数精度数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS复化积分数值积分公式与多项式插值有很大的关系。因此Runge现象的存在，使得我们不能用太多的积分点计算。采用与插值时候类似，我们采用分段、低阶的方法数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)( 12)()(2)(311iiixxfhxfxfhxfii1031

19、11031)( 12)(21)()(21 )( 12)()(2)(niiniiniiiinfhbfxfafhfhxfxfhfT误差误差做等距节点，niihaxnabhi, 0,复化梯形公式数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS由均值定理知)( )( .,.,102nfftsbabaCfnii)( 12)()( )(12)( 12)(2323fnabfabhfnhfEn可以看出，复化梯形公式是收敛的。如果节点不等距，还可以做复化积分吗？怎么处理？数 学 系University of S

20、cience and Technology of ChinaDEPARTMENT OF MATHEMATICS)(2880)2()()(4)(62)()4(522122222iiiixxfhxfxfxfhxfii10)4(5112101210)4(522122)(2880)2()()(2)(4)(3 )(2880)2()()()(62)(miimiimiiniiiiinfhbfxfxfafhfhxfxfxfhfS误差误差做等距节点，mnniihaxnabhi2;, 0,复化Simpson公式数 学 系University of Science and Technology of ChinaDE

21、PARTMENT OF MATHEMATICS由均值定理知)(180)()(2880)()(2880)2()()4(45)4(45)4(5fnabfmabfmhfEn可以看出，复化Simpson公式是收敛的。数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS定义定义 若一个积分公式的误差满足若一个积分公式的误差满足 且且C 0，则，则称该公式是称该公式是 p 阶收敛阶收敛的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例：例：计算计算dxx 10142 解：解： )1

22、()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfxffSkk8kxk 其中其中= 3.141592502运算量基运算量基本相同本相同数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSLab03 复化积分1.分别编写用复化Simpson积分公式和复化梯形积 分公式计算积分的通用程序2.用如上程序计算积分51( )sin( )I fx dx取节点xi , i=0,N，N 为 2k,k=0,1,12 ，并计算误差

23、，同时给出误差阶3.简单分析你得到的数据数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS误差阶：记步长为h时的误差为e,步长为h/k时的误差为ek则，相应的误差阶为：)ln()ln(keeokk数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSSample Output ( represents a space)复化梯形积分，误差和误差阶为k=0,0.244934066848e00k=1,0.5

24、34607244904 , 1.90.复化Simpson积分，误差和误差阶为k=1,0.244934066848e00k=2,0.534607244904e-01 , 4.01.数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS 函数变化有急有缓，为了照顾变化剧烈部分的误差，我们需要加密格点。对于变化缓慢的部分，加密格点会造成计算的浪费。以此我们介绍一种算法，可以自动在变化剧烈的地方加密格点计算，而变化缓慢的地方，则取稀疏的格点。积分的自适应计算积分的自适应计算数 学 系University o

25、f Science and Technology of ChinaDEPARTMENT OF MATHEMATICS先看看事后误差估计事后误差估计以复化梯形公式为例)( 12)()()(2fhabfTfIn)( 212)()()(22fhabfTfInn等分区间2n等分区间近似有：)( )( ff)()(31)()(22fTfTfTfInnn)()(151)()(22fSfSfSfInnn类似，复化Simpson公式数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS自适应计算记为复化一次，2次的

26、Simpson公式,21baSbaS0控制dxxffIba)()(求数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS,21baSbaS15,21baSbaS,)(2baSfI是215,2,22152,2,2121bbaSbbaSbaaSbaaS数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)()(31)()(22fTfTfTfInnn由前面的事后误差估计式，)()()(31)()(222fS

27、fTfTfTfInnnn则，这启发我们，可以用低阶的公式组合后成为一个高阶的公式。)()()(151)()(222fCfSfSfSfInnnn类似，Romberg积分积分数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS记为以步长为h的某数值积分公式，有)(hI)()()(mmhochhIfImmhohchIfI22)2()(12)(22)(mhIhIhIfI 1222)(mhIhIhIfI数 学 系University of Science and Technology of ChinaDEP

28、ARTMENT OF MATHEMATICS有如下的Euler-Maclaurin定理)(12)()2()2()2(222)()()()1(mmmmmmhOhIhIhIhI若)()()(2)(mmhOhIfI为2m阶公式，则Romberg 积分就是不断地用如上定理组合低阶公 式为高阶公式，进而计算积分 Romberg 算法：算法： ? ? ? T1 =)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R1 =)0(3T S2 =)1(1T C1 =)0(2T C2 =)1(2T S4 =)2(1T数 学 系University of Science a

29、nd Technology of ChinaDEPARTMENT OF MATHEMATICS重积分的计算 badcdxdyyxf),( 在微积分中，二重积分的计算是用化为累次积分的方法进行的。计算二重数值积分也同样采用累次积分的计算过程。简化起见，我们仅讨论矩形区域上的二重积分。对非矩形区域的积分，大多可以变化为矩形区域上的累次积分。a,b,c,d 为常数，f 在D 上连续。将它变为化累次积分 dcbabadcbadcdydxyxfdxdyyxfdxdyyxf),(),(),(首先来看看复化梯形公式的二重推广数 学 系University of Science and Technology

30、of ChinaDEPARTMENT OF MATHEMATICS做等距节点，x轴，y轴分别有：ncdkmabh,dcdyyxf),(先计算),(21),(),(21),(110nnjjdcyxfyxfyxfkdyyxf，将x作为常数，有再将y作为常数，在x方向，计算上式的每一项的积分),(21),(),(212),(210110000yxfyxfyxfhdxyxfmmiiba),(21),(),(212),(21110nmmininbanyxfyxfyxfhdxyxf二重积分的复化梯形公式二重积分的复化梯形公式数 学 系University of Science and Technology

31、 of ChinaDEPARTMENT OF MATHEMATICS 1111110111101111),(),(21),(21 ),(21),(),(21 ),(),(njmijinjjmjnjjmmijijnjbajbanjjyxfhyxfyxfhyxfyxfyxfhdxyxfdxyxf数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS 1111,1111111101111000000),( ),( ),(),(),(),(21 ),(),(),(),(41),(njmijijinjmiji

32、njjmnjjminimiimmnbadcyxfchkyxfyxfyxfyxfyxfyxfyxfyxfyxfhkdxdyyxf系数，在积分区域的四个角点为1/4，4个边界为1/2，内部节点为1),(),(12)()(222222fykfxhabcdfE误差数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS类似前面有： minjjijibadcyxfhkdxdyyxf00,),(),(记TnTmvvvVuuuU31,34,32,34,32,34,31,31,34,32,34,32,34,31,10

33、10jijivu ,二重积分的复化二重积分的复化Simpson公式公式做等距节点，x轴，y轴分别有：ncdkmabh,m，n为偶数数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS),(),(180)()(444444fykfxhabcdfE误差数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSGauss型积分公式 Newton-Cotes 积分公式，可以知道n为偶数时，n+1个点数值积分公式有n

34、+1阶精度。是否有更高的代数精度呢？n个点的数值积分公式，最高可以到多少代数精度？本节会解决这个问题。数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS例：在两点数值积分公式中，如果积分点也作为未知量，则有4个未知量，可以列出4个方程： （以f(x)在-1,1为例）032021111133113001122112001111001110dxxxaxadxxxaxaxdxxaxadxaa可解出：31,31, 1, 11010 xxaa数值积分公式)31()31(11fffdx具有3阶代数精度，比梯

35、形公式1阶代数精度高数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSbaiiniiinbadxxlaxfafIdxxffI)(, )()()()(0证明：取)()()()()(222120 xxxxxxxxpnn易知：0)(0)(xpIxpIn也就是说，数值积分公式，对一个2n+2阶的多项式是有误差的，所以，n+1个点的数值积分公式不超过2n+1阶n 个积分点的数值积分公式，最高2n1阶定理如何构造如何构造最高阶精度的公式？数 学 系University of Science and Tech

36、nology of ChinaDEPARTMENT OF MATHEMATICS一般性，考虑积分：0)(,)()()(xWdxxfxWfIba称为权函数定义两个可积函数的内积为：badxxgxfxWgf)()()(),(两个函数正交，就是指这两个函数的内积为0数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS利用Schmidt 正交化过程，nxx, 1变为正交基0010( )( )( ),( )( )( )( )( ),( )nninniiiigxfxfx g xgxfxg xg x g x)(

37、,),(),(10 xpxpxpn就可以将多项式基函数数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS以n阶正交多项式的n个零点为积分点的数值积分公式有2n1阶的代数精度Gauss点Gauss积分，记为Gn(f)证明：bannndxxWxxxxxffIfIfE)()(,)()()(21,21xxxxfn若 f 为 2n1 次多项式，则为 n1 次多项式又，)(),(xpxnn仅差一个常数（零点相同）0)(fE1)(nnPfp具有一个很好的性质：nfIfGn),()(数 学 系Universit

38、y of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS(2)求出pn(x)的n n个零点x1 , x2 , xn 即为Gauss点. (1)求出区间a,b上权函数为W(x)的正交多项式pn(x) .(3)计算积分系数 Gauss型求积公式的构造方法数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)()(),()(,()(00001xpxpxpxpxxxp)()(),()(,()()(),()(,()(111120000222x

39、pxpxpxpxxpxpxpxpxxxp5321141151121142xxdxxdxxdxxdxxx解解 按 Schemite 正交化过程作出正交多项式: 的2点Gauss公式.求积分dxxfx)(112例：0( )1p x 数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS故两点Gauss公式为 积分系数为31)(11212211121dxxxxxxdxxlxA112212211211( )3xxAx lx dxxdxxx)()()(535331112ffdxxfxP2(x)的两个零点为 ,

40、532531xx数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS 区间-1,1上权函数W(x)=1的Gauss型求积公式,称为Gauss-Legendre求积公式求积公式,其Gauss点为Legendre多项式的零点. (1) Gauss-Legendre求积公式求积公式公式的Gauss点和求积系数可在数学用表中查到 .几种几种Gauss型求积公式型求积公式由因此,a,b上权函数W(x)=1的Gauss型求积公式为batabbaxdttabbafabdxxf)2)()()22(2)(11ban

41、iiixabbafAabdxxf1)22(2)(数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSnxkAknxkAk10260.93246951420.66120938650.2386191861036076157300.467913934620.5773502692130.774596669200.55555555560.888888888970.94910791230.74153118560.405845151400.12948496620.27970539150

42、.38183005050.417959183740.86113631160.33998104360.34785484510.652145154980.96028985650.79666647740.5255324099010122853630.22238103450.31370664590.362683783450.90617984590.538469310100.23692688510.47862867050.5688888889数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHE

43、MATICS 区间0,)上权函数W(x)=e-x的Gauss型求积公式,称为Gauss-Laguerre求积公式求积公式,其Gauss点为Laguerre多项式的零点. (2) Gauss-Laguerre求积公式求积公式公式的Gauss点和求积系数可在数学用表中查到 .由所以,对0, + )上权函数W(x)=1的积分,也可以构造类似的Gauss-Laguerre求积公式:00)()(dxxfeedxxfxx01)()(niixixfeAdxxfi数 学 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSnxkAknxkAk20.58588643763.41421356230.85355339050.146446609450.26356031971.41340305913.5964