微分方程数值求解

1、数值积分引言引言引言引言引言n由定积分定义iniiniiniiiiiinbainiixxfsSxxxxfsbxxxaxfdxxf0011000)()3()()2(.)1()(lim)(求和近似分割引言baniiixnxinidxxfxfSxx)()(limlimmax)4(0001求极限的近似值。也是加权和，是权系数，其中由此想到机械求积公式bainiiiiniiinnbadxxfxfxfAAfRxfAfRxfAxfAxfAdxxf)()()()()(.)()()(0011005.1 Newton-Cotes求积公式),()(12)()(2)()(21)()()( 3 bafabfRfRbfa

2、fabdxbxaxfdxbfabaxafbabxdxxfTTbababa）（其中所以n由Lagrange插值,任何一的函数 都可以近似的表示成其中)(xfy )()()(xRxLxfnn插值多项式。的是LagragexfyxlxLnjjjn)()()(0n为简便起见,取节点为等分n现在关键是求njjhaxnabhj,.,2 , 1 , 0,dxxRdxxLdxxfbanbaban)()()( banjiiijibajnjjnjnjjnjbajjnjbajbanjjjbandxxxxxabdxxlabCxfCabxfdxxlababydxxldxyxldxxL0)(0)(0001)(1)()()

3、()(1)()()()(因此就归结为求权 nnjiinjiinnjiibanjiiijinjijiijdtitijnhdtijitnhdxxxxxabCntbxtaxhijxxhitxxhdtdxthaxihaxnjjhaxnabh000000)()(1111;0,)(,)(,.,2, 1 ,0,因此。时时，知由21211)0()!11 (! 11) 1(212) 1(11) 1()!01 (! 01) 1(1)()!( !) 1(1021011)1(11021001)1(000 tdttCtdttCndtitjnnjnnjiijn时，仅有两个节点：当616461)2(21)2(3141)2(

4、)2(41)2)(1()!02(!02)1(2)2(2)2(120232022002)2(0CCttdtttdtttCn，同理可得时当n以此类推得Cotes系数表:Newton Cotes积分公式)()()()()()()()(,)(),.,2, 1 ,0(,)(0)(0fRfCabdxxRdxxLdxxfxRxLxflagragexfxfxfnjkhaxnabhnbabaxfnjjnjbababannnnnjjjj则称插值多项式，即的为节点作以记等分，取区间等分上的连续函数，将是设定义的截断误差。为，（其中称系数；是其中积分公式。 banbanbannnnjiijnnjdxxfxxxxxxx

5、dxxxxxxfdxxfnfRCotesdtitjnnjCCotesNewton)().()()()(,.,)()()!1(1)()!( !)1(1010)1(00)(常用的几个积分公式n梯形公式(n=1),()(12)()()(2)(21)(21)()(,21 31110bafabfRbfafabbfafabfTfRfTdxxfCCTbaT且。则因为nSimpson公式(n=2),()(2880)()()2(4)(6)(61,64,61)4(5)2(2)2(1)2(0bafabfRbfbafafabfSfRfSdxxfCCCSbaS且所以因为nNewton公式(n=3)。其中且所以，因为3)

6、()2(3)(3)(8)(8183,83,81)3(3)3(2)3(1)3(0abhbfhafhafafabfNfRfNdxxfCCCCbaNnCotes公式(n=4)。其中且所以，因为4)(7)2(32)2(12)(32)(7(90)(907,90329012,9032,907)4(4)4(3)4(2)4(1)4(0abhbfhafhafhafafabfCfRfCdxxfCCCCCbaC例题756931.075693.0)21351334131(81;4469.0)21231411(61;57 .011212169314718.02ln112121ICotesINewtonISimpsonI

7、dxxILeibnizNewtondxxI公式得由公式由公式由）（由梯形公式公式得解：由。计算5.1.2 Newton-Cotes公式截断误差及代数精度公式截断误差及代数精度,().(2)(1()!1()().(2)(1)(2()!2()(,)(2880)(,)(1 .5.10)1(20)2(3)4(54banabhndtntttnfhndtntttntnfhfRCotesNewtonnbafabfRSimpsonbaCxfnnnnnnnS其中，为奇即数）为偶数）公式的截断误差为一般地对任意的截断误差为公式则设定理,)(483840)(4)6(7bafabfRCotesnC求积公式的截断误差得

8、令几个常用的求积公式的代数精度1.T 公式的代数精度度公式具有一次的代数精所以时当时当TdxxfbaabbfafabfTabxdxxdxxfxxfdxxfbaabbfafabfTabxxdxdxxfxxfbabababababababa)()(2)()(2)(3131)()()()(2)()(2)(2121)()(22333222222. S公式的代数精度成立所以时当)()(21)24(6)()2(4)(6)(21)()(2222fSdxxfabbbaaabbfabfafabfSabxdxdxxfxxfbababa精确成立即时当)()(31)242(6)2(4(6)()2(4)(6)(31)(

9、)(33222223322fSdxxfabbabaabbbaaabbfabfafabfSabdxxdxxfxxfbababa精确成立即时当)()(41)(236)33(21(6)2(4(6)()2(4)(6)(41)()(4432233322333334423fSdxxfabbabbaaabbbabbaaaabbbaaabbfabfafabfSabdxxdxxfxxfbababan因此S-公式具有三次代数精度。n同理可得N-公式具有三次代数精度，C-公式具有五次代数精度。0).(1(,2).(1()(2222202dumuuuhfRmtumnndtnttthdxxfRmmmnban则为偶数，所

10、以令因为5.2 复化求积公式n将区间a,b适当分割成若干个字区间，对每个子区间使用求积公式，构成所谓的复化求积公式，这是提高积分精度的一个常用的方法。5.2.1定步长复化求积公式n1.复化梯形求积公式 bbaa2)()2(2)(4)()2(22)2()(22)2()()(2)(bfbafafabbfbafbabafafbahTbfafbahTn一般地将a,b区间n等分，则)()(2),.2 , 1(,),.2 , 1 , 0(,211jjjjjjxfxfhSTnjxxnjjhaxbah公式有使用对每个子区间n所以)() )(2)()(2) )(2)()(2)()()()(.)()()()(2)

11、 )()(2,)(1112111111hTjhafbfafhxfbfafhbfbfbfxfxfxfxfafhxfxfhShfRSdxxfnnjnjjnnjjjnjjTbanjj而n而)(12)()(12,)(1)(,)(),()(12, 2 31 11 3fhabnfhhfRfnfbabaxfxxfhhfRTnjjjjjnjjT故）使得（则必存在一点区间上的二阶连续函数是若定步长复化梯形求积公式算法ShSSnjjhafSSbfafSnabhhTnban输出计算输入.3)3(;,.2, 1)()2(;2)()(;)1()(.2;,.1n2.复化Simpson公式类似于梯形公式：,)()()(4)

12、(6,.2 , 1,12111hfRSdxxfxfxfxfhSnjxxSnjjbajjjjjj则上有在每个子区间njnjjjnjnjjjnnnnnnjjjjnjjxfxfbfafhxfxfbfafhfffffffffffhxfxfxfhS112111212112231121012111)()(22)()(3)(2)(4)()(64.446)()(4)(6（而),()(2880)(),()(2880,)()21(2)(2)()(3) )()21(22)()(3)4(411)4(5111bafhabxxfhhfRhShjafjhafbfafhjhafhjafbfafhjjinjiSnnjnjnj定

13、步长复化Simpson求积公式算法SShSnjjhafhjafSSbfafSnabhhSnban输出计算输入. 3;3)3();,.2 , 1()()21(2)2(;2)()(;) 1 ()(. 2;,. 1例题至少取多少？试问划分数要使得截断误差不超过公式计算用复化如果用复化梯形公式和对于定积分例nxdxI5201021,Simpsonsin1 . 2 . 5254102148sin224,)(sin)2(1202)(12)(,52 2 2nnnhfRnfhabhfRTT解得即解：由截断误差有。取故解得即由6, 1 . 531. 02,31. 01021sin22880,)(sin)(288

14、002)(2880)(,54)4(4)4(4nnhhhfRhfhabhfRSS5.2.2变步长求积公式变步长梯形求积公式算法5.2.1变步长梯形求积算法);(21)2();)21() 1 ()(. 3;2)()(; 1)(. 2;,. 112121nnjnnnHTThjafhHhTbfafhTabhnhTMba计算计算输入3;2;2. 6thenif . 5;thenif . 421212返回计算；输出划分数过大，停止输出TThhnnMnTTT变步长Simpson求积公式)2(2)2(31)2()(2)(31)21(2)(2)()(31)21(2)(2)()(3)(222111hHhThShH

15、hThjafhjhafbfafhhjafjhafbfafhhSSimpsonnnnnnnjnjnjn所以公式由复化重复上述过程。将分点加密一倍，定积分的近似值，否则为时，取当程序实现的基本思想：)2()()2()2()2()2()()()(22222hShShShShHhThShHhTnnnnnnnnn变步长Sinmpson求积算法);(31);)21()2(;2;2);(21) 1 ()2(),2(),2(. 3);2(31);2()2(;2)()(; 1) 1 ()(),(),(. 2;,. 122112222111nnjnnnnnnnnnnHTShjafhHhhnnHTThShHhTHT

16、SbahfHbfafhTabhnhShHhTMba计算计算输入。返回计算；输出分划数过大，停止停止计算；输出3;. 6thenif. 5,thenif. 42121212SSTTMnSSS14)()2(4)(14)()2(4)(2222hShShCCoteshThThSSimpsonnnnnnn公式同理也可以推出复化有如下关系公式与复化梯形公式复化例题。于公式计算，使其误差小用变步长计算，并估计误差；的复化用公式；用计算51010)3(5)2() 1 (1SimpsonSimpsonnSimpsonxdxI所以有：则节点）（)5 , 4 , 3 , 2 , 1 , 0(02 . 051)2(6

17、9444. 03625)2123141 (61)1 ()21(4)0(611iihxnabhfffIi69315. 0)218 . 116 . 114 . 112 . 11(221(301)1118 . 0116 . 0114 . 0112 . 011(2)111011(5161)(5hSI5445)4(103333. 1120)2 . 0(242880, 1 , 024124)(hhfRxxxfS所以因为5.3 Romberg求积公式n5.3.1外推法基本思想 以较小的计算量为代价,达到提高数值结果的精度是外推法的中心思想.展式起截断误差为由。逼近用中心差商设Taylorxfhhxfhxfh

18、GhxhxCxfk)()2()2()(,2,2)(22)2()2(.)2()2()2()(2)()()(),()()2(),.,2 , 1()!12(2)()(.)()(224422222)12(2224422hEhhhhGxfhhhOhGxfhOhETaylarhxfhkjjxfhEhhhhGxfkkkjjjkk代得用上式中。且展式得的证明可以从无关的常数。是与其中)()()(),.3 , 2() 14(31,14)()(4)(,14)()2(4)()()()(.)2()()()(4412)1(2111111226644hOhGxfkjrhEhEhEhGhGhGhGxfhEhrhrhrhGxfhGxfjjjkk此时。其中,.2 , 1 , 014)2()(4)()()(),()()(14)2()(4)()()2()(11106221122211mhGhGhGhGhGhOhGxfhGhGhGhGhGhGmmmmm依此类推有：这时构造和同理可以有外推法。上述过程称为而RichardsonhOhGxfmm),()()(225.3.2 Romberg求积算法求积算法Romberg求积算法 5.4 Gauss求积公式5.4.1 正交多项式正交多项式