《欧拉公式的改进ppt课件.ppt》由会员分享，可在线阅读，更多相关《欧拉公式的改进ppt课件.ppt（22页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、1 Eulers Method 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy+ + )1,., 0(),(111 = =+ += =+ + + +niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公
2、式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(+=iiiyxyR)()(322hOxyih+ + = =即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 Hey! Isnt the leading term of the local truncation error of Eulers method ? Seems that we can make a good use of it )(22ihxy 1 Eulers Method 梯形公式梯形公式 /* trapezoid formula
3、 */ 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 = =+ + += =+ + + +niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii=+ 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxy
4、xyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy+ + 1,., 1),(211 = =+ += = + +niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy= = = )()(311hOyxyRiii= = = =+ + +需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-st
5、ep method */。方方 法法 1 Eulers Method显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, let
6、s make it possible. 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+ += =+ +Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+ +iy),(),(2111+ + + + + += =iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到
7、它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。 )1,., 0(),(,),(211 = =+ + + += =+ + +niyxfhyxfyxfhyyiiiiiiii1 Eulers Method2 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到
8、沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii+ + += = = + + += =+ +斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的
9、前提假设下，使得 )(iixyy = =)()(311hOyxyRiii= = = =+ + +Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii+ + + += =+ + += =)()()(2hOxyphxyii+ + + + = =将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+ + += = =+ + += =+ + ),(),(),(),(),(),()(yxfyxfyx
10、fdxdyyxfyxfyxfdxdxyyxyx+ += =+ += = = Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii+ + + + + + += =+ + + + + + + += =+ + 2 Runge-Kutta MethodStep 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii+ + + + + + += =+ + )()(2)()()(321hO
11、xyhxyhxyxyiiii+ + + + + += =+ +要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii=+21,1221= = =+ +p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121= = = = p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij
12、 ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 + + + + + + += =+ + + += =+ + += = =+ + + + += =mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Met
13、hod */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii+ + += =+ + += =+ + += = =+ + + + += =+ +2 Runge-Kutta Method注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数
14、)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。HW: p.202 #1, 22 Runge-Kutta MethodLab 15. Runge-Kutta (Order Four) Use Runge-Kutta method of order four to approximate the solution of the i
15、nitial-value problem, , and (1)You are supposed to write a functionvoid RK4 ( double (*f)( ), double a, double b, double y0, int n, FILE *outfile)to approximate the solution of Problem (1) with y= f, x in a, b, and the initial value of y being y0. Output the approximating values of y on the n+1 equa
16、lly spaced grid points from a to b to outfile.InputThere is no input file. Instead, you must hand in your function in a *.h file. The rule of naming the *.h file is the same as that of naming the *.c or *.cpp files.Output ( represents a space) For each test case, you are supposed to print n+1 lines,
17、 and each line is in the following format: fprintf( outfile, “%8.4f%12.8en”, x, y );),(yxfy = = ,bax 0)(yay= =2 Runge-Kutta MethodSample Judge ProgramSample Judge Program#include #include #include 98115001_15.h double f0(double x, double y) return (y x*x+1.0); void main( ) FILE *outfile = fopen(out.
18、txt, w);int n;double a, b, y0; a = 0.0; b = 2.0; y0 = 0.50; n = 10;RK4(f0, a, b, y0, n, outfile);fprintf(outfile, n);fclose(outfile);Sample Output Sample Output ( ( represents a space) represents a space) 0.00005.00000000e 0010.20008.29293333e 0010.40001.21407621e+0000.60001.64892202e+0000.80002.127
19、20268e+0001.00002.64082269e+0001.20003.17989417e+0001.40003.73234007e+0001.60004.28340950e+0001.80004.81508569e+0002.00005.30536300e+0003 收敛性与稳定性收敛性与稳定性 /* Convergency and Stability */ 收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算
20、法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 = = = 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)(= =欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 + += =+ += =+ +0)1 (yhyii + += =对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0+ += =+ += =ehhh=+ /10)1(lim)(0ixxyeyi= = 3 Convergency and Stability 稳定性稳定性 /* Stabili
21、ty */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 = = = = 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧欧拉隐式拉隐式欧拉显式欧拉显式 节点节点 xixey30= 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.0000
22、2.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!3 Convergency and Stability定义定义若某算法在计算过程中任一步产
23、生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */yy = = 常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定
24、区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。000yy = = h h= =h3 Convergency and Stability例：例：考察显式欧拉法考察显式欧拉法011)1(yhyhyyiiii+ + + += =+ += = 000yy = = 011)1(yhyii+ + + += =01111)1( + + + + + += = = =iiiihyy由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh = =1|1| + + h0-1-2
25、ReImg例：例：考察隐式欧拉法考察隐式欧拉法11+ + + += =iiiyhyy iiyhy = =+ +11101111 + + + = =iih可见绝对稳定区域为：可见绝对稳定区域为：1|1| h210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。法的好。3 Convergency and Stability例：例：隐式龙格隐式龙格-库塔法库塔法 = =+ + + + += =+ + + += =+ +),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方
26、法的绝对稳定阶方法的绝对稳定区域为区域为 + + += =+ += =+ +)2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定HW: p.202 #64 线性多步法线性多步法 /* Multistep Method */用用若干若干节点处的节点处的 y 及及 y 值的值的线性组合线性组合来近似来近似y(xi+1)。).(.110111101kikiiikikiiiffffhyyyy + + + + + + + + + + + += = 其通式可写为：其通式
27、可写为：),(jjjyxff = =当当 1 0 时，为时，为隐式公隐式公式式; 1=0 则为则为显式公式显式公式。 基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxfy =,1+iixx + += = + +1)(,()()(1iixxiidxxyxfxyxy只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通，则可通过过 近似近似y(xi+1) 。而。而选用不同近似式选用不同近似式 Ik，可得到不，可得到不同的计算公式同的计算公式。 + + 1)(,(iixxkdxxyxfIkiiIyy+ += =+ +14 Multistep Method 亚
28、当姆斯显式公式亚当姆斯显式公式 /* Adams explicit formulae */利用利用k+1 个节点上的被积函数值个节点上的被积函数值 构造构造 k 阶牛顿阶牛顿后插后插多项式多项式 , 有有kiiifff ,.,11, 0, )( + +thtxNik + + + += =+ +1010)()()(,(1dthhtxRdthhtxNdxxyxfikxxikiiNewton插值余项插值余项dthtxNhyyikii)(101+ + += = + +/* 显式计算公式显式计算公式 */局部截断误差为：局部截断误差为： + += = = =+ + +1011)()(dthtxRhyxy
29、Rikiii例：例：k=1 时有时有)()(11 + += = + += =+ +iiiiiifftfftfhtxN + + + += = + + += =10111)3(2)(iiiiiiiiffhydtfftfhyydththtdxyfdhRxxi)1(!21)(,(1022+ += = )(1253iyh = =4 Multistep Method注：注：一般有一般有 ，其中，其中Bk 与与yi+1 计算公式计算公式中中 fi , , fi k 各项的各项的系数系数均可查表得到均可查表得到 。 )()2(2ikkkiyhBR + + += =10123k21231223245521121
30、62459125243724912583720251fifi 1fi 2fi 3BkMisprint on p.204常用的是常用的是 k = 3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式)9375955(243211 + + + + + += =iiiiiiffffhyy4 Multistep Method 亚当姆斯隐式公式亚当姆斯隐式公式 /* Adams implicit formulae */利用利用k+1 个节点上的被积函数值个节点上的被积函数值 fi+1 , fi , , fi k+1 构造构造 k 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项
31、式完全类似地可得到一系列隐式公式隐式公式，并有，并有 ，其中，其中 与与 fi+1 , fi , , fi k+1 的系数亦可查表得到。的系数亦可查表得到。)()2(2ikkkiyhBR + + += =kB10123k21 21125249211282419121 245 241121 241 72019 fi+1fifi 1fi 2Bk常用的是常用的是 k = 3 的的4阶亚当姆斯隐式公式阶亚当姆斯隐式公式)5199(242111 + + + + + + += =iiiiiiffffhyy小于小于Bk较同阶显较同阶显式式稳定稳定4 Multistep Method 亚当姆斯预测亚当姆斯预测
32、-校正系统校正系统 /* Adams predictor-corrector system */Step 1: 用用Runge-Kutta 法法计算前计算前 k 个初值；个初值；Step 2: 用用Adams 显式显式计算计算预测预测值；值；Step 3: 用同阶用同阶Adams 隐式隐式计算计算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必须相同。式的精度必须相同。通常用通常用经典经典Runge-Kutta 法法配合配合4阶阶Adams 公式公式。 Hey! Look at the local truncation error of the explicit and implic
33、it Adams methods: and Dont you think theres something you can do?)(720251)5(5iyh )(72019)5(5iyh 4阶阶Adams隐式公式的截断误差为隐式公式的截断误差为)(72019)()5(511iiiyhyxy = = + + +4阶阶Adams显式公式的截断误差为显式公式的截断误差为)(720251)()5(511iiiyhyxy = = + + +当当 h 充分小充分小时，可近似认为时，可近似认为 i i ，则：，则：19251)()(1111 + + + + +iiiiyxyyxy)(270251)(11
34、11+ + + + + + + iiiiyyyxy)(27019)(1111+ + + + + iiiiyyyxyPredicted value pi+1Modified value mi+1Corrected value ci+1Modified final value yi+1外推技术外推技术 /* extrapolation */4 Multistep Method Adams 4th-Order predictor-corrector AlgorithmTo approximate the the solution of the initial-value problemAt (N+1
35、) equally spaced numbers in the interval a, b.Input: endpoints a, b; integer N; initial value y0 .Output: approximation y at the (N+1) values of x.Step 1 Set h = (b a) / N ; x0 = a; y0 = y0; Output ( x0, y0 );Step 2 For i = 1, 2, 3 Compute yi using classical Runge-Kutta method; Output ( xi , yi );St
36、ep 3 For i = 4, , N do steps 4-10Step 5 ; /* predict */ Step 6 ; /* modify */Step 7 ; /* correct */Step 8 ; /* modify the final value */Step 9 Output ( xi+1 , yi+1 ); Step 10 For j = 0, 1, 2, 3 Set xi = xi+1 ; yi = yi+1 ; /* Prepare for next iteration */Step 11 STOP. 0)(,),(yaybxayxfy= = = =24/ )9375955(3211 + + + + + += =iiiiiiffffhyp270/ )(25111iiiipcpm + += =+ + +24/ )519),(9(21111 + + + + + + + += =iiiiiiifffmxfhyc270/ )(191111+ + + + + = =iiiipccy应为应为( ci+1 pi+1 ), 但因但因ci+1 尚未尚未算出，只好用算出，只好用( ci pi )取代之。取代之。