计算机数值方法 第五章 常微分方程数值解法lz

第五章常微分方程数值解法(NumericalMethodsforOrdinaryDifferentialEquations)引言引言考虑一阶常微分方程的初值问题(Initial-ValueProblem):0)(],[),(yaybaxyxfdxdy只要f(x,y)在x∈[a,b]上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意定义在[a,b]上的y1(x)和y2(x)都成立，则上述常微分方程存在唯一解。|||),(),(|2121yyLyxfyxf引言数值解法就是要计算出解函数y(x)在一系列节点a=x0x1…xn=b处的近似值y0y1…yn节点间距称为步长，通常采用等距节点，即取hi=h(常数)。)1,...,0(1nixxhiii解决：数值解法的一个基本特点是“步进式”，即求解时顺着节点排列的次序一步步地向前推进。0)(],[),(yaybaxyxfdxdy单步：yk-1yk多步：yk-p…yk-2，yk-1yk注意：与“迭代法”区别第五章常微分方程数值解法欧拉方法（Euler’sMethod）显式Euler公式121()()(()2)jjjjy'xyxyxyhh''由两点公式求导数，在[xj,xj+1]子区间上有：其中j[xj,xj+1]12()())()(2jjjjyxyxhy'hy''x2()(,()2)jjjjyxhfxyhy''代入方程有(,)dyfxydx1(,)jjjjyyhfxy显式Euler公式x0x1y(x1)y(x0)h显式Euler公式的误差1(,)jjjjyyhfxy21112()()22()jjjjjeyxhhy''y''xy局部截断误差显式Euler公式定义在假设yj=y(xj)，即第i步计算是精确的前提下，考虑的截断误差ej+1=y(xj+1)yj+1称为局部截断误差/*localtruncationerror*/。12()()(,()2)jjjjjyxyxhfxhy''y（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题]1,0[1)0(2xyyxyy解：用显式Euler法求解，有：)2(),(y1iiiiiiiiyxyhyyxhfy…依次下去计算结果见P1690100012111220y()10.1(1)1.11220.1y()1.10.1(1.1)1.19181.1xyhyyxyhyy隐式Euler公式121(1()(2)())jjjjy'xyxyxhyh''由两点公式求导数，在[xj,xj+1]子区间上有：其中j[xj,xj+1]112()()()()2jjjjyxyhxhy'xy''211()(()2,)jjjjyxhfxyhy''代入方程有(,)dyfxydx111(,)jjjjyyhfxy221111()()2)2(jjjjjhhy'e'xxyyy''局部截断误差隐式Euler公式是一个关于yj+1的方程,要从中解出yj+1x0x1（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题]1,0[1)0(2xyyxyy计算过程见P170用隐式Euler法求解，有：)2(),(y111111iiiiiiiiyxyhyyxhfy解：从中解出yi+1，有：21i18(1)y2(1)iiiyyhhxh显隐结合的预测校正系统——避免求解方程（predictor-correctormethod）Step1:先用显式欧拉公式作预测，算出预测值Step2:再用隐式欧拉公式作校正，得到校正值)1,...,0(),(,11niyxfhyxhfyyiiiiii111(,)jjjjyyhxyf1(,)jjjjyhfxyy写成一个公式为：y0＝y(a)－y1－y1－y2－y2…－yn－yn计算顺序：（P170）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题]1,0[1)0(2xyyxyy……计算过程见P170解：用预测校正系统求解，有：i11i1112y()2y()iiiiiiiixyhyyxyhyy01000110112y()1.12y()1.0918xyhyyxyhyyy0=112111221222y()1.18272y()1.1763xyhyyxyhyyp阶精度若某算法的局部截断误差e(h)满足：e(h)=O(hp+1)，即有：e(h)/hp+1=c(常数），则称该算法有p阶精度。定义欧拉法的局部截断误差：欧拉法具有1阶精度。2111()2()jjjjeyhy''xy2111()2()jjjjeyxhy''y显式：隐式：梯形公式（trapezoidformula）)1,...,0()],(),([2111niyxfyxfhyyiiiiii注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式。)()(311hOyxyRiii从显式、隐式Euler法和Euler局部截断误差来看似乎可以有如下式子——梯形公式：1121(2,())jjjjjyyhfxyhy''12(,)2()jjjjjyyhfxyhy''隐式Euler显式Euler基于数值积分的求解思想111()()()(,)kkkkxxkkxxyxyxyxdxfxydx11()()(,)kkxkkxyxyxfxydx求出积分的近似表达，也就求出了y(xk)1(,)kkxxfxydxhxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)基于数值积分的求解思想用矩形求积公式代替11()()(,)kkxkkxyxyxfxydxhxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)11(,)kkhfxy111(,)kkkkyyhfxy1(,)kkkkyyhfxy(,)kkhfxy左矩形显式Euler公式右矩形隐式Euler公式基于数值积分的求解思想用梯形求积公式代替11()()(,)kkxkkxyxyxfxydxhxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)11(,)(,)2kkkkhfxyfxy111(,)(,)2kkkkkkhyyfxyfxy注：梯形公式是隐式公式显式、隐式Euler的平均梯形公式具有2阶精度33()()12khehf''Oh余项方法显式欧拉隐式欧拉梯形公式计算简单精度低稳定性好精度低,计算不便精度提高计算不便Can’tyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?Letmetry！Euler公式及梯形公式总结改进的Euler法——预测校正系统（modifiedEuler’smethod）注：可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的稳定性高于显式欧拉法。误差比较：梯形法≈改进Euler法显、隐Euler法Step1:先用显式欧拉公式作预测，算出预测值Step2:再用梯形公式作校正，得到校正值111(,)kkkkfyyyhx111(,)(,)2kkkkkkhyyfxyyfx举例(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题]5.0,0[1)0(1xyxyy解：计算过程参见P173Euler法111(1)(1)2kkkkkkhyyyxyx11(2)()22kkkkhyhxxhyh111(1)kkkkyyhyx梯形法解出yk得：举例(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题]5.0,0[1)0(1xyxyy解：计算过程参见P173改进Euler法111(1)(1)2kkkkkkhyyyxxy111(1)kkkkyhyxySimpson公式（Simpsonformula）注：Simpson公式是一个隐式公式Simpson公式是多步法：yk-1,yk-2=yk5(5)5()()()90hehyOh应用simpson求积公式：22()()(,)kkxkkxyxyxfxydx221(4)6kkkkkhyyfffSimpson公式——4阶精度P198习题6用Euler法和改进的Euler法解方程，h＝0.1注意改题目为0=x=0.4第五章作业欧拉方法（Euler’sMethod）第五章常微分方程数值解法龙格－库塔法（Runge－KuttaMethod）二阶Runge－Kutta法——改进的Euler法1(,)iiiiyhfxyy111(,)(,)2iiiiiihyyfxyfxy考察改进的欧拉法，可以将其改写为：K1K21iyhKxi+h1121211122(,)(,)iiiiiiyyhKKKfxyKfxhyhK步长一定是一个h吗？斜率一定取K1K2的平均值吗？()iyx龙格-库塔法（Runge－KuttaMethod）首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得)(iixyy)()(311hOyxyRiiiStep1:将K2在(xi,yi)点作Taylor展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii)()()(2hOxyphxyii将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx龙格-库塔法（Runge－KuttaMethod）Step2:将K1K2代入第1式，得到)()()()()]()()([)(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiiiStep3:将yi+1与y(xi+1)在xi点的泰勒展开作比较)()()()(322211hOxyphxyhyyiiii)()(2)()()(321hOxyhxyhxyxyiiii要求，则必须有：)()(311hOyxyRiii龙格-库塔法（Runge－KuttaMethod）21,1221p这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意:就是改进的欧拉法,也就是2阶R-K法。21,121p:为获得更高的精度，应该如何进一步推广？三阶龙格-库塔法三、四阶R-K法思想：利用各种显化方法，对以上公式显化11))(,)((iixxii