南大数值分析课件第九章常微分方程数值解.ppt

《南大数值分析课件第九章常微分方程数值解.ppt》由会员分享，可在线阅读，更多相关《南大数值分析课件第九章常微分方程数值解.ppt（18页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第九章第九章常微分方程数值解常微分方程数值解/*Numerical Methods for Ordinary Differential Equations*/考虑考虑一阶一阶常常微分方程的微分方程的初值问题初值问题/*Initial-Value Problem*/:只要只要f(x,y)在在a,b R1上连续，且关于上连续，且关于y满足满足Lipschitz条条件件，即存在与，即存在与x,y 无关的常数无关的常数L使使对任意定义在对任意定义在a,b上的上的y1(x)和和y2(x)都成立，则上述都成立，则上述IVP存在唯一解存在唯一解。要计算出解函数要计算出解函数y(x)在一系列节点在一系列节点a

2、=x0 x1xn=b 处的近似值处的近似值节点间距节点间距为步长，通常采用为步长，通常采用等距节点等距节点，即取，即取hi=h(常数常数)。1欧拉方法欧拉方法/*Eulers Method*/欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数记为记为在在假假设设yi=y(xi)，即即第第 i步步计计算算是是精精确确的的前前提提下下，考考虑虑的的截截断断误误差差Ri=y(xi+1)yi+1称称为为局局部部截截断断误误差差/*local truncation error*/。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p阶精度。阶精度。欧拉法

3、的局部截断误差：欧拉法的局部截断误差：欧拉法具有欧拉法具有1阶精度。阶精度。Ri的的主项主项/*leading term*/亦亦称为称为欧拉折线法欧拉折线法/*Eulers polygonal arc method*/1EulersMethod 欧拉公式的改进：欧拉公式的改进：隐式欧拉法隐式欧拉法/*implicit Euler method*/向后差商近似导数向后差商近似导数x0 x1)(,()(1101xyxfhyxy+)1,.,0(),(111=+=+niyxfhyyiiii由于未知数由于未知数yi+1同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐

4、式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显式显式/*explicit*/欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即即隐式欧拉公式具有隐式欧拉公式具有1阶精度。阶精度。Hey!Isnttheleading termofthelocal truncation errorofEulersmethod?Seemsthatwecanmakeagooduseofit1EulersMethod梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法的显、

5、隐式两种算法的平均平均注：注：的确有局部截断误差的确有局部截断误差，即梯形公式具有即梯形公式具有2阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。中点欧拉公式中点欧拉公式/*midpoint formula*/中心差商近似导数中心差商近似导数x0 x2x1假设假设，则可以导出，则可以导出即中点公式具有即中点公式具有2阶精度。阶精度。需要需要2个初值个初值y0和和y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双

6、步法/*double-step method*/，而前面的三种算法都是而前面的三种算法都是单步法单步法/*single-step method*/。方方法法 1EulersMethod显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高,显式显式多多一个初值一个初值,可能影响精度可能影响精度Cantyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?Doyouthinkitposs

7、ible?Well,callmegreedyOK,letsmakeitpossible.改进欧拉法改进欧拉法/*modified Eulers method*/Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法/*predictor-corrector method*/。可以证明该算法具有可以证明该算法具有2阶精度，同时可以看到它是个阶精度，

8、同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。1EulersMethod2龙格龙格-库塔法库塔法/*Runge-Kutta Method*/建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。考察改进的欧拉法，可以将其改写为：考察改

9、进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2的的平均值平均值吗？吗？步长一定是一个步长一定是一个h吗吗？2Runge-KuttaMethod首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有使得到的算法格式有2阶阶精度，即在精度，即在的前提假设下，使得的前提假设下，使得 Step 1:将将K2在在(xi,yi)点作点作Taylor展开展开将将改进欧拉法推广为：改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将K2代入第代入第1式，得到式，得到2Runge-KuttaMethodStep 3:

10、将将yi+1与与y(xi+1)在在xi 点点的的泰勒泰勒展开作比较展开作比较要求要求，则必须有：，则必须有：这里有这里有个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。2Runge-Ku

11、ttaMethod).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：2Runge-KuttaMethod注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算Ki的值，即计算的值，即计算f的的值。值。Butcher于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：

12、度阶数的关系：753可可达到的最高精度达到的最高精度642每步须算每步须算Ki 的的个数个数由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h取小取小。HW:p.202#1,22Runge-KuttaMethodLab15.Runge-Kutta(OrderFour)UseRunge-Kuttamethodoforderfourtoapproximatethesolutionoftheinitial-valueprob

13、lem,and(1)YouaresupposedtowriteafunctionvoidRK4(double(*f)(),doublea,doubleb,doubley0,intn,FILE*outfile)toapproximatethesolutionofProblem(1)withy=f,xina,b,andtheinitialvalueofybeingy0.Outputtheapproximatingvaluesofyonthen+1equallyspacedgridpointsfromatobtooutfile.InputThereisnoinputfile.Instead,youm

14、usthandinyourfunctionina*.hfile.Theruleofnamingthe*.hfileisthesameasthatofnamingthe*.cor*.cppfiles.Output(representsaspace)Foreachtestcase,youaresupposedtoprintn+1lines,andeachlineisinthefollowingformat:fprintf(outfile,“%8.4f%12.8en”,x,y);2Runge-KuttaMethodSample Judge ProgramSample Judge Program#in

15、clude#include#include98115001_15.hdoublef0(doublex,doubley)return(y x*x+1.0);voidmain()FILE*outfile=fopen(out.txt,w);intn;doublea,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.000000

16、00e 0010.20008.29293333e 0010.40001.21407621e+0000.60001.64892202e+0000.80002.12720268e+0001.00002.64082269e+0001.20003.17989417e+0001.40003.73234007e+0001.60004.28340950e+0001.80004.81508569e+0002.00005.30536300e+0003收敛性与稳定性收敛性与稳定性 /*Convergency and Stability*/收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固

17、固定定的的x=xi=x0+i h，当当h0(同时同时i )时有时有yiy(xi)，则称该算法是则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为欧拉欧拉公式为公式为对对任意固定的任意固定的x=xi=i h，有，有 3ConvergencyandStability 稳定性稳定性/*Stability*/例：例：考察初值问题考察初值问题在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.4

18、0.5精确解精确解改进欧拉法改进欧拉法 欧欧拉拉隐式隐式欧拉欧拉显式显式 节点节点xi1.0000 2.00004.0000 8.00001.6000 101 3.2000 1011.00002.5000 10 16.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7Whatiswrong?!AnEngineercomplains:Maththe

19、oremsaresounstablethatasmallperturbationontheconditionswillcauseacrashontheconclusions!3ConvergencyandStability若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定的的/*absolutely stable*/。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程/*test equation*/常数，可以常数，可以是复数是复数当当步长取为步长取为h时，将某

20、算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差，则若此误差以后逐步衰减，就称该算法，则若此误差以后逐步衰减，就称该算法相对于相对于绝对稳定绝对稳定，的全体构成的全体构成绝对稳定区域绝对稳定区域。我们。我们称称算法算法A比算法比算法B稳定稳定，就是指，就是指A的绝对稳定区域比的绝对稳定区域比B的的大大。h h=h3ConvergencyandStability例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。3ConvergencyandStability例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式14阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中其中2阶方法阶方法的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定HW:p.308#5