Euler法与改进Euler法

1、常微分方程数值解法-欧拉法、改进欧拉法和四阶龙格库塔法常微分方程数值解法常微分方程数值解法 常微分方程主要有常微分方程主要有: : (1) (1)变量可分离的方程变量可分离的方程 (2)(2)一阶线性微分方程一阶线性微分方程( (贝努利方程贝努利方程) ) (3) (3)可降阶的一类高阶方程可降阶的一类高阶方程 (4)(4)二阶常系数齐次微分方程二阶常系数齐次微分方程 (5)(5)二阶常系数非齐次微分方程二阶常系数非齐次微分方程 (6)(6)全微分方程全微分方程常微分方程数值解法常微分方程数值解法l主要内容：l一、 引引 言言l二、二、 建立数值解法的常用方法建立数值解法的常用方法l三、三、

2、EulerEuler方法方法l四、四、 几何意义几何意义l五、五、 EulerEuler方法的误差估计方法的误差估计l六、六、 改进欧拉法改进欧拉法l七、四阶龙格库塔法七、四阶龙格库塔法l七、程序设计要求七、程序设计要求主要内容主要内容 许多实际问题的数学模型是微分方程或微分许多实际问题的数学模型是微分方程或微分方程的定解问题方程的定解问题, ,如物体运动如物体运动, ,电路震荡电路震荡, ,化学反化学反映及生物群体的变化等映及生物群体的变化等. . 能用解析方法求出精确解的微分方程为数不能用解析方法求出精确解的微分方程为数不多多, ,而且有的方程即使有解析解而且有的方程即使有解析解, ,也可

3、能由于解的也可能由于解的表达式非常复杂而不易计算表达式非常复杂而不易计算, ,因此有必要研究微因此有必要研究微分方程的数值解法分方程的数值解法 一、引一、引 言言重点 研究一阶常微分方程的初值问题的数值解)1()(),(:00bxayxyyxfdxdy 其一般形式为其一般形式为假定:.),(),(:,),(论知论知这样由常微分方程的理这样由常微分方程的理）条件）条件满足利普希茨（满足利普希茨（且关于且关于连续连续函数函数yyLyxfyxfLipschitzyyxf .)(存在并且唯一存在并且唯一解解xyy )1(初值问题初值问题常微分方程数值解法常微分方程数值解法bxxxxxaxyNn 210

4、)(在一系列离散节点在一系列离散节点求解求解所谓数值解法，就是寻所谓数值解法，就是寻初值问题数值解的提法为定数，为定数，定定如不特别说明，总是假如不特别说明，总是假称为步长。称为步长。相邻两个节点的间距相邻两个节点的间距上的近似值上的近似值)2 , 1 , 0(,1210ihhxxhyyyyyinnnNn.2 , 1 , 0,0nnhxxn常微分方程数值解法常微分方程数值解法建立微分方程数值解法,首先要将微分方程离散化.一般采用以下几种方法:(1) 用差商近似导数11:() ,()nnnnyy xyy x进一步令dxdyhyynn1)(,11,nnnnnnxyxfxxxyxyyxdxdynn

5、二、建立数值解法的常用方法二、建立数值解法的常用方法(2) (2) 用数值积分近似积分用数值积分近似积分 1)(,()()(1nnxxnndxxyxfxyxy即即)(,)(:11nnnnxyyxyy 令令进一步进一步),()(,(11nnxxnnyxfhdxxyxfyynn 实际上是矩形法实际上是矩形法宽宽高高), 1 , 0(),(11 ndxyxfdxdxdynnnnxxxx常用方法常用方法(3) (3) 用用TaylorTaylor多项式近似并可估计误差多项式近似并可估计误差)( ! 2)( )()( ! 2)( )()()(221nnnnnnnxyhxhyxyyhxhyxyhxyxy

6、)(,)(:11nnnnxyyxyy 令令进一步进一步),(1nnnnyxhfyy )( max2)(211xyhyxybxann 常用方法常用方法)1()(),(:00bxayxyyxfdxdy 式为式为已知初值问题的一般形已知初值问题的一般形用差商近似导数用差商近似导数dxdyhyynn 1问题转化为问题转化为,.)3 , 2 , 1 , 0()(),(001 nxyyyxhfyynnnnEulerEuler方法的迭代公式方法的迭代公式 三、三、EulerEuler方法方法111:(,)(0)1nnnnyyhKKf xyy作等价变换,有,.)3 , 2 , 1 , 0()(),(001 n

7、xyyyxhfyynnnn1(,)nnKf xy令令EulerEuler方法方法 已知，必有切线方程。已知，必有切线方程。及及由于由于斜率斜率的切线（存在！），则的切线（存在！），则出发取解曲线出发取解曲线由由 , , , , 000000,0000yxyxfyxfyxxy yyxdxdy ),()( 0000,0000yxfxxyyxxxyydxdy ：由点斜式写出切线方程由点斜式写出切线方程 四、几何意义四、几何意义）（：，可由切线算出，可由切线算出，则，则为为等步长等步长0001101, yxhfyyyhxxh 2 1 0 , )(n1n1n，）（：点的值点的值在在逐步计算出逐步计算出

8、nyxhfyyxxyynn注意：这是“折线法”而非“切线法”除第一个点是曲线切线外，其他点不是！Y=y(x)ab1x2x几何意义几何意义五、Euler方法的误差估计为简化分析，先考虑计算一步所产生的误差，即假设为简化分析，先考虑计算一步所产生的误差，即假设)(nnxyy 是精确的，是精确的， 估计误差估计误差111)(nnnyxyR这种误差称为这种误差称为局部截断误差局部截断误差。 估计截断误差的主要方法是估计截断误差的主要方法是Taylor展开法，即将函数展开法，即将函数)(xy在在nx处展开：处展开： )(2)()()(2nnnnxyhxyhxyhxy取一次取一次Taylor多项式近似函数

9、，得多项式近似函数，得 )()(1hxyxynn)(2)()(2yhxyhxynn )(21)(,()(2yhxyxhfxynnn 12 )(21),( nnnnnxxyhyxhfy得得Euler方法的方法的局部截断误差公式局部截断误差公式为为 )(21)(2111yhyxyRnnn 结论：结论：上式说明上式说明Euler公式公式的的局部截断误差局部截断误差为为)(2hO它的精度很差。它的精度很差。 一般很少用它来求近似值，但是一般很少用它来求近似值，但是Euler法法却体现了数值方法的基本思想。却体现了数值方法的基本思想。 定义定义8.1 如果某种数值方法的局部截断误差为如果某种数值方法的局

10、部截断误差为1()pO h，则称该方法是，则称该方法是p阶方法或阶方法或具有具有p阶精度阶精度。显然。显然p越大，方法的越大，方法的精度越高。精度越高。 注注：Euler方法具有一阶精度，因此它的精度不高。方法具有一阶精度，因此它的精度不高。 六六改进的改进的EulerEuler方法方法改进的改进的Euler方法方法)1()(),(:00bxayxyyxfdxdy 式为式为已知初值问题的一般形已知初值问题的一般形利用数值积分将微分方程离散化利用数值积分将微分方程离散化得梯形公式得梯形公式: :),(),(211 nnnnyxfyxfh 1)(,()()(1nnxxnndxxyxfxyxy),(

11、),(2111 nnnnnnyxfyxfhyy解决方法：有的可化为显格式，但有的不行解决方法：有的可化为显格式，但有的不行梯形方法为隐式算法梯形方法为隐式算法改进的改进的Euler方法方法 ,kyxfyxfhyyyxhfyynEulerknnnnnknnnnn210 ,2, 210)(11)1(1)0(1，对，对法结合，形成迭代算法法结合，形成迭代算法与与梯形公式比欧拉法精度高一些梯形公式比欧拉法精度高一些, ,但计算量较大但计算量较大 实际计算中只迭代一次，这样建立的预测实际计算中只迭代一次，这样建立的预测校正系统称作校正系统称作改进的欧拉公式改进的欧拉公式。改进的改进的Euler方法方法)

12、,(),(2111 nnnnnnyxfyxfhyy 00121211)(),(),()(2:yxyhKyhxfKyxfKKKhyynnnnnn作等价变换作等价变换),(),(2),(1111 kkkkkkkkkkyxfyxfhyyyxhfyy校正校正预测预测改进的改进的Euler方法方法二、改进的Euler法 梯形方法虽然提高了精度，但算法复杂，计算梯形方法虽然提高了精度，但算法复杂，计算量大。如果实际计算时精度要求不太高，用梯形公量大。如果实际计算时精度要求不太高，用梯形公式求解时，每步可以迭代一次，由此导出一种新的式求解时，每步可以迭代一次，由此导出一种新的方法方法改进改进Euler法法。

13、这种方法实际上就是。这种方法实际上就是将将Euler公式与梯形公式结合使用公式与梯形公式结合使用：先用：先用Euler公式公式求求 1ny的一个初步近似值的一个初步近似值 1ny，称为，称为预测值预测值，预测值，预测值 的精度可能很差，再用的精度可能很差，再用梯形公式校正梯形公式校正求得近似值求得近似值 1ny ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyy即即 校正预测改进改进Euler法法 亦称为由亦称为由Euler公式公式和和梯形公式梯形公式得到的得到的 预测校正预测校正(Predictor-Corrector)系统。系统。 )(21),( ),(1qp

14、npnnqnnnpyyyyhxhfyyyxhfyy为便于上机编程，常改写成为便于上机编程，常改写成l改进改进EulerEuler方法计算框图方法计算框图开始开始00, ,xy h b输输入入1n 1000001112(,)(,)()pcpcpxxhyyhf xyyyhf xyyyy11,xy输输出出1xb结结束束01011,nnxxyyY YN N 1)0()10(2)1 . 0(yxyxyyh步步长长求求解解初初值值问问题题例例解解yxyyxf/2),( （1）用）用Euler方法方法得算式为得算式为 1)0()2(1yyxyhyynnnnn（2）用）用改进的改进的Euler方法方法得算式为

15、得算式为 )2(0.1nnnnpyxyyy )1 . 0(2( 1 . 0pnpnqyxyyy )(211qpnyyy 七、七、 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的基本思想是从单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一点出发，以某一斜率沿直线达到斜率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变点。欧拉法及其各种变形所能达到的最高精度为形所能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(21

16、21121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的平均值吗？的平均值吗？步长一定是一个步长一定是一个h 吗？吗？首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将

17、改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhx

18、yhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在无穷多个解。所有满足上式的格式统称为存在无穷多个解。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式。塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均

19、为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶经典龙格阶经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 2 Runge-Kutta Metho

20、d注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长采用低阶算法而将步长h h 取小。取小。lRunge-KuttaRunge-