第二章數值微分和數值積分

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)無關，與積分區(qū)間和積分點有關數 學 系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 積分若節(jié)點可以自由選取，則，一個自然的辦法就是取等距節(jié)點。對區(qū)間做等距分割。該數值積分稱為Newton-Cotes積分數 學 系University of Science and Technology of ChinaDEPARTMENT OF MATHE

14、MATICSniihaxnabhi,0,dtntititttnininhdthinintititttdxxlaninninthaxbaii00)() 1)(1() 1()!( !) 1( ) 1()!( !)() 1)(1() 1()()()(niiCaba設節(jié)點步長)( 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誤差誤差做等距節(jié)點，niihaxnabhi, 0,復化梯形公式數 學 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS由均值定理知)( )( .,.,102nfftsbabaCfnii)( 12)()( )(12)( 12)(2323fnabfabhfnhfEn可以看出，復化梯形公式是收斂的。如果節(jié)點不等距，還可以做復化積分嗎？怎么處理？數 學 系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誤差誤差做等距節(jié)點，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取節(jié)點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等分區(qū)間2n等分區(qū)間近似有：)( )( 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則，這啟發(fā)我們，可以用低階的公式組合后成為一個高階的公式。)()()(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),( 在微積分中，二重積分的計算是用化為累次積分的方法進行的。計算二重數值積分也同樣采用累次積分的計算過程。簡化起見，我們僅討論矩形區(qū)域上的二重積分。對非矩形區(qū)域的積分，大多可以變化為矩形區(qū)域上的累次積分。a,b,c,d 為常數，f 在D 上連續(xù)。將它變?yōu)榛鄞畏e分 dcbabadcbadcdydxyxfdxdyyxfdxdyyxf),(),(),(首先來看看復化梯形公式的二重推廣數 學 系University of Science and Technology

30、of ChinaDEPARTMENT OF MATHEMATICS做等距節(jié)點，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系數，在積分區(qū)域的四個角點為1/4，4個邊界為1/2，內部節(jié)點為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公式公式做等距節(jié)點，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個點的數值積分公式，最高可以到多少代數精度？本節(jié)會解決這個問題。數 學 系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變?yōu)檎换?010( )( )( ),( )( )( )( )( ),( )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)求出區(qū)間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 區(qū)間-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 區(qū)間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