当前位置:首页 > 办公文档 > 会议纪要 > 计算方法5.常微分方程的数值解法概要
2020/1/1815.3龙格-库塔法/*Runge-KuttaMethod*/考察改进的欧拉法,可以将其改写为:),(),(2121121211hKyhxfKyxfKKKhyyiiiiii斜率一定取K1K2的平均值吗?步长一定是一个h吗?单步递推法的基本思想是从(xi,yi)点出发,以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。建立高精度的单步递推格式。2020/1/182首先希望能确定系数1、2、p,使得到的算法格式有2阶精度,即在的前提假设下,使得)(iixyy)()(311hOyxyRiiiStep1:将K2在(xi,yi)点作Taylor展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii)()()(2hOxyphxyii将改进欧拉法推广为:),(),(][12122111phKyphxfKyxfKKKhyyiiiiii),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyxStep2:将K2代入第1式,得到)()()()()]()()([)(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii2020/1/183Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较)()()()(322211hOxyphxyhyyiiii)()(2)()()(321hOxyhxyhxyxyiiii要求,则必须有:)()(311hOyxyRiii21,1221p这里有个未知数,个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。21,121p注意到,就是改进的欧拉法。Q:为获得更高的精度,应该如何进一步推广?2020/1/184构造高阶单步法的直接方法•由Taylor公式:•当h充分小时,略去Taylor公式余项,并以yi、yi+1分别代替y(xi)、y(xi+1),得到差分方程:•其局部截断误差为:121()(1)()()()'()()...()()2!!(1)!iippppiiiiyxyxhhhhyxhyxyxyxypp),(!...),('!2),()1(21iippiiiiiiyxfphyxfhyxhfyy)()!1(~)()1(111ppiiyphyxy差分方程为p阶方式,上述方式称为Taylor方式。注:利用Tayler公式构造,不实用,高阶导数f(i)不易计算。2020/1/185RungeKutta方法1.基本思想因为其中K=f(,y())称为y(x)在[xi,xi+1]上的平均斜率。若取K1=f(xi,y(xi))——Euler公式取K2=f(xi+1,y(xi+1))——向后Euler公式一阶精度取——梯形公式二阶精度猜想:若能多预测几个点的斜率,再取其加权平均作为K,可望得到较高精度的数值解,从而避免求f的高阶导数。11()()(,())()(,())()iixiixiiyxyxfxyxdxyxhfyyxhK)(2121KK2020/1/1862.RK公式其中Kj为y=y(x)在xi+ajh(0aj1)处的斜率预测值。aj,bjs,cj为特定常数。111112),,(),(jssjsijijiipjjjiipjKbhyhaxfKyxfKKchyy2020/1/1873.常数的确定确定的原则是使精度尽可能高。以二阶为例:(希望y(xi+1)–yi+1=O(hp)的阶数p尽可能高)首先:),(),()(12122122111hKbyhaxfKyxfKKcKchyyiiiiii)()(!21)(')()(321hOxyhxhyxyxyiiii2020/1/188另一方面:将K2在(xi,yi)处展开。K2=f(xi,yi)+a2hfx'(xi,yi)+b21hK1fy'(xi,yi)+O(h2)可得:yi+1=yi+hc1f(xi,yi)+hc2f(xi,yi)+h2c2[a2fx'(xi,yi)+b21K1fy'(xi,yi)]+O(h3)=yi+h(c1+c2)f(xi,yi)+c2a2h2[fx'(xi,yi)+(b21/a2)f(xi,yi)fy'(xi,yi)]+O(h3)(希望)...)),()(),(),((000000yxfyyxxyxfyyxxf)()()(')(322221hOxyhacxycchyiii2020/1/189希望:ei+1=y(xi+1)–yi+1=O(h3).则应:特例:a2=1c1=c2=1/2,b21=1,得2阶R-K公式:改进欧拉公式。12112212221abaccc),(),()(2121211hKyhxfKyxfKKKhyyiiiiii2020/1/1810希望:ei+1=y(xi+1)–yi+1=O(h3).则应:特例:c1=0c2=1,a2=1/2,b21=1/2,得:称为中点公式。12112212221abaccc)2,2(),(12121KhyhxfKyxfKhKyyiiiiii2020/1/18114.最常用的R-K公式——标准4阶R-K公式算法:),()2,2()2,2(),()22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiiii输入a,b,n,y0h=(b-a)n,x0=afori=1,i=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0,y02020/1/1812Matlab代码:functionRunge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;fori=1:nK1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy1=y0+h*(K1+2*K2+2*K3+K4)/6;y0=y1end;end;functionf=f(x,y)f=x+y;end;输入a,b,n,y0h=(b-a)n,x0=afori=1,i=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0,y02020/1/1813例:用标准4阶R-K公式求:的数值解。取h=0.2,并与标准解y=2ex–x–1比较。解:因为f(x,y)=x+y,从由得:1)0(10)()(yxxyxxy3423121432112.02.022.022.022.022.0)22(62.0KyxKKyxKKyxKyxKKKKKyyiiiiiiiiii2020/1/1814xi+1K1K2K3K4yi+1y(xi+1)0.211.21.221.4441.2428001.2428060.41.4428001.6870801.7115081.9851021.5836361.5836490.61.9836362.2820002.3118362.6460032.0442132.0442380.82.6442133.0086343.0450763.4532282.6510422.6510821.03.4510423.8961463.9406574.4391733.4365032020/1/1815注:龙格-库塔法的主要运算在于计算Ki的值,即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系:753可达到的最高精度642每步须算Ki的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n由于龙格-库塔法的导出基于泰勒展开,故精度主要受解函数的光滑性影响。对于光滑性不太好的解,最好采用低阶算法而将步长h取小。深入研究龙格-库塔法请看此处!2020/1/18161.Runge-Kutta法的一般形式2.2阶Runge-Kutta方法3.经典Runge-Kutta方法4.R-K-Fehhlberg方法5.隐式R-K方法6.变步长方法龙格-库塔法深入研究2020/1/18171.Runge-Kutta法的一般形式Runge-Kutta法是用区间上若干个点上的导数的线性组合得到平均斜率,以提高方法的阶。其一般形式为:11Lnniiiyyhk1222133331132211,,1,2,,(,)(,)(,())nnnnnniininiijjjkfxchychakiLkfxykfxchychkkfxchychakak2020/1/1818式(9.11)称L级p阶Runge-Kutta方法(简称R-K法)。当L=1就是欧拉法,当L=2时为改进的欧拉法。1111,1,1Liiiijijca。111()()LnnniiiTyxyxhk其中它的局部截断误差是(9.11)2.2级2阶Runge-Kutta方法令L=2,则11122()nnyyhkk2020/1/18203.经典Runge-Kutta方法我们可以构造出一族3级3阶,一族4级4阶和一族5级4阶等R-K方法。最常用的4级4阶是如下的经典R-K方法:112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyykkkkkfxyhhkfxykhhkfxykkfxhyhk2020/1/18214.R-K-Fehhlberg方法R-K-Fehhlberg方法是在R-K方法的基础上引进误差和步长控制的办法。即利用5阶R-K法112456166656285619213512825564305055iiyyKKKKK估计4阶R-K的局部误差,其中5431151410421972565104821625KKKKyyii)329323,83()41,4(),(213121KKyhxhfKKyhxhfKyxhfKiiiiii2020/1/1822)401141041859256553442278,2()410485451336808216439,()219772962197720021971932,1312(543216432153214KKKKKyhxhfKKKKKyhxhfKKKKyhxhfKiiiiii注:R-K-Fehhlberg比4阶R-K方法具有更大的优越性,他是计算稳定,高精度的方法,他的显著优点是,每一步仅需计算f的6个值;若用4阶R-K-L与5阶R-K-L在一起使用,则每步需要计算f的10个值。推荐使用!2020/1/18235.隐式R-K方法类似于显式R-K公式,稍加改变,就得到隐式
本文标题:计算方法5.常微分方程的数值解法概要
链接地址:https://www.777doc.com/doc-3150831 .html