Euler法与改进Euler法

常微分方程数值解法-欧拉法、改进欧拉法和四阶龙格库塔法常微分方程数值解法常微分方程主要有:(1)变量可分离的方程(2)一阶线性微分方程(贝努利方程)(3)可降阶的一类高阶方程(4)二阶常系数齐次微分方程(5)二阶常系数非齐次微分方程(6)全微分方程常微分方程数值解法主要内容：一、引言二、建立数值解法的常用方法三、Euler方法四、几何意义五、Euler方法的误差估计六、改进欧拉法七、四阶龙格库塔法七、程序设计要求主要内容许多实际问题的数学模型是微分方程或微分方程的定解问题,如物体运动,电路震荡,化学反映及生物群体的变化等.能用解析方法求出精确解的微分方程为数不多,而且有的方程即使有解析解,也可能由于解的表达式非常复杂而不易计算,因此有必要研究微分方程的数值解法一、引言重点研究一阶常微分方程的初值问题的数值解)1()(),(:00bxayxyyxfdxdy其一般形式为假定:.),(),(:,),(论知这样由常微分方程的理）条件满足利普希茨（且关于连续函数yyLyxfyxfLipschitzyyxf.)(存在并且唯一解xyy)1(初值问题常微分方程数值解法bxxxxxaxyNn210)(在一系列离散节点求解所谓数值解法，就是寻初值问题数值解的提法为定数，定如不特别说明，总是假称为步长。相邻两个节点的间距上的近似值)2,1,0(,,,,1210ihhxxhyyyyyinnnNn.2,1,0,0nnhxxn常微分方程数值解法建立微分方程数值解法,首先要将微分方程离散化.一般采用以下几种方法:(1)用差商近似导数11:(),()nnnnyyxyyx进一步令dxdyhyynn1)(,11,nnnnnnxyxfxxxyxyyxdxdynn二、建立数值解法的常用方法(2)用数值积分近似积分1))(,()()(1nnxxnndxxyxfxyxy即)(,)(:11nnnnxyyxyy令进一步),())(,(11nnxxnnyxfhdxxyxfyynn实际上是矩形法宽高),1,0(),(11ndxyxfdxdxdynnnnxxxx常用方法(3)用Taylor多项式近似并可估计误差)(''!2)(')()(''!2)(')()()(221nnnnnnnxyhxhyxyyhxhyxyhxyxy)(,)(:11nnnnxyyxyy令进一步),(1nnnnyxhfyy)(''max2)(211xyhyxybxann常用方法)1()(),(:00bxayxyyxfdxdy式为已知初值问题的一般形用差商近似导数dxdyhyynn1问题转化为,...)3,2,1,0()(),(001nxyyyxhfyynnnnEuler方法的迭代公式三、Euler方法111:(,)(0)1nnnnyyhKKfxyy作等价变换,有,...)3,2,1,0()(),(001nxyyyxhfyynnnn1(,)nnKfxy令Euler方法已知，必有切线方程。及由于斜率的切线（存在！），则出发取解曲线由,,,,000000,0000yxyxfyxfyxxyyyxdxdy),()(0000,0000yxfxxyyxxxyydxdy：由点斜式写出切线方程四、几何意义）（：，可由切线算出，则为等步长0001101,yxhfyyyhxxh210,)(n1n1n，，，）（：点的值在逐步计算出nyxhfyyxxyynn注意：这是“折线法”而非“切线法”除第一个点是曲线切线外，其他点不是！yχ0χ1χ2χ3χY=y(x)ab1x2x几何意义五、Euler方法的误差估计为简化分析，先考虑计算一步所产生的误差，即假设)(nnxyy是精确的，估计误差111)(nnnyxyR这种误差称为局部截断误差。估计截断误差的主要方法是Taylor展开法，即将函数)(xy在nx处展开：)(2)()()(2nnnnxyhxyhxyhxy取一次Taylor多项式近似函数，得)()(1hxyxynn)(2)()(2yhxyhxynn)(21))(,()(2yhxyxhfxynnn12)(21),(nnnnnxxyhyxhfy得Euler方法的局部截断误差公式为)(21)(2111yhyxyRnnn结论：上式说明Euler公式的局部截断误差为)(2hO它的精度很差。一般很少用它来求近似值，但是Euler法却体现了数值方法的基本思想。定义8.1如果某种数值方法的局部截断误差为1()pOh，则称该方法是p阶方法或具有p阶精度。显然p越大，方法的精度越高。注：Euler方法具有一阶精度，因此它的精度不高。六改进的Euler方法改进的Euler方法)1()(),(:00bxayxyyxfdxdy式为已知初值问题的一般形利用数值积分将微分方程离散化得梯形公式:)],(),([211nnnnyxfyxfh1))(,()()(1nnxxnndxxyxfxyxy)],(),([2111nnnnnnyxfyxfhyy解决方法：有的可化为显格式，但有的不行梯形方法为隐式算法改进的Euler方法,,,kyxfyxfhyyyxhfyynEulerknnnnnknnnnn210,,2,210)(11)1(1)0(1，，，，对法结合，形成迭代算法与梯形公式比欧拉法精度高一些,但计算量较大实际计算中只迭代一次，这样建立的预测—校正系统称作改进的欧拉公式。改进的Euler方法)],(),(2111nnnnnnyxfyxfhyy00121211)(),(),()(2:yxyhKyhxfKyxfKKKhyynnnnnn作等价变换))~,(),((2),(~1111kkkkkkkkkkyxfyxfhyyyxhfyy校正预测改进的Euler方法二、改进的Euler法梯形方法虽然提高了精度，但算法复杂，计算量大。如果实际计算时精度要求不太高，用梯形公式求解时，每步可以迭代一次，由此导出一种新的方法——改进Euler法。这种方法实际上就是将Euler公式与梯形公式结合使用：先用Euler公式求1ny的一个初步近似值1ny，称为预测值，预测值的精度可能很差，再用梯形公式校正求得近似值1ny)],(),([2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy即校正预测改进Euler法亦称为由Euler公式和梯形公式得到的预测－校正(Predictor-Corrector)系统。)(21),(),(1qpnpnnqnnnpyyyyhxhfyyyxhfyy为便于上机编程，常改写成改进Euler方法计算框图开始00,,,xyhb输入1n1000001112(,)(,)()pcpcpxxhyyhfxyyyhfxyyyy11,xy输出1xb结束01011,nnxxyyYN1)0()10(2)1.0(yxyxyyh步长求解初值问题例解yxyyxf/2),(（1）用Euler方法得算式为1)0()2(1yyxyhyynnnnn（2）用改进的Euler方法得算式为)2(0.1nnnnpyxyyy))1.0(2(1.0pnpnqyxyyy)(211qpnyyy七、龙格-库塔法/*Runge-KuttaMethod*/建立高精度的单步递推格式。单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii斜率一定取K1K2的平均值吗？步长一定是一个h吗？首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得)(iixyy)()(311hOyxyRiiiStep1:将K2在(xi,yi)点作Taylor展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii)()()(2hOxyphxyii将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyxStep2:将K2代入第1式，得到)()()()()]()()([)(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiiiStep3:将yi+1与y(xi+1)在xi点的泰勒展开作比较)()()()(322211hOxyphxyhyyiiii)()(2)()()(321hOxyhxyhxyxyiiii要求，则必须有：)()(311hOyxyRiii21,1221p这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。21,121p注意到，就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？其中i(i=1,…,m)，i(i=2,…,m)和ij(i=2,…,m;j=1,…,i1)均为待定系数，确定这些系数的步骤与前面相似。)...,(......),(),(),(]...[1122112321313312122122111mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii§2Runge-KuttaMethod注：龙格-库塔法的主要运算在于计算Ki的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：753可达到的最高精度642每步须算Ki的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。Runge-Kutta方法的设计思想设法在[xn,xn+1]区间内多预报几个点的斜率值，利用这些斜率值，将他们加权平均作为平均斜率的近似，有可能构造出更高精度的计算格式分别用欧拉法、改进欧拉法、四阶龙格库塔法及Dsolve方法求解并对比.八、程序设计要求1)0()10(2)1.0(yxyxyyh步长求解初值问题11)0()10(32)1.0(2yxyxyh步长求解初值问题2解常微分方程的初值问题)10(2.1)0(2xyxydxdy用欧拉法求解分别取步长h=0.1和h=0.05。3