数值计算课件——第六章常微分方程初值问题的数值解法n

1、常微分方程初值问题常微分方程初值问题的数值解法的数值解法第第6章章引言引言在实际问题中，常需要求解微分方程在实际问题中，常需要求解微分方程(如发电机运动方程如发电机运动方程)。只有简单的和典型的微分方程可以求出解析解，而在只有简单的和典型的微分方程可以求出解析解，而在实际问题中的微分方程往往无法求出解析解。实际问题中的微分方程往往无法求出解析解。在高等数学中我们见过以下常微分方程：在高等数学中我们见过以下常微分方程： 0)(),(yaybxayxfy-(1) )(,)(),(0ayyaybxayyxfy-(2)一阶常微分方程一阶常微分方程 nybyyaybxayyxfy)(,)(),(0-(3

2、)(1),(2)式称为式称为初值问题初值问题,(,(3)式称为式称为边值问题边值问题 2002212210012111)(),()(),(yxyyyxfyyxyyyxfy-(4)另外另外, ,在实际应用中还经常需要求解常微分方程组在实际应用中还经常需要求解常微分方程组: 本课程主要研究问题本课程主要研究问题一阶常微分方程一阶常微分方程(1)的数值解法的数值解法,我们首先介绍初值问题我们首先介绍初值问题(1)(1)的解存在的条件的解存在的条件定理定理 只要只要 f (x, y) 连续，且关于连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与即存在与 x, y 无关的常数无关的常数

3、 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则初值问都成立，则初值问题（题（1）存在唯一解存在唯一解。| ),(),(|2121yyLyxfyxf bxxxxan 210), 2 , 1()(nkyxykk 的的近近似似值值上上函函数数值值的数值解的数值解就是问题就是问题而而)1(), 2 , 1(nkyk 上上的的一一系系列列离离散散点点在在区区间间就就是是求求未未知知函函数数,)(baxy（通常采用（通常采用等距节点）等距节点）对于问题对于问题(1)(1) 0)(),(yaybxayxfy要求它的要求它的数值解数值解常微分方程数值解公式的推

4、导常微分方程数值解公式的推导 求初值问题数值解的方法是求初值问题数值解的方法是步进法步进法，即从已知的初值，即从已知的初值y0出发，通过出发，通过一定的计算一定的计算求求y1 ，然后由，然后由y1或或y0和和y1求出求出y2 ，依次计算到依次计算到yn ，即在计算出，即在计算出yk后计算后计算yk+1 ，这时有，这时有单步法单步法：计算：计算yk+1时，只利用时，只利用yk多步法多步法：计算：计算yk+1时，用到时，用到yk, yk-1, yk-2,常微分方程数值解公式常微分方程数值解公式的主要推导方法的主要推导方法泰勒展开泰勒展开利用差商利用差商利用数值积分法利用数值积分法1、泰勒展开的求解

5、思路、泰勒展开的求解思路：将将 按泰勒级数展开按泰勒级数展开 hxyxykk 1 kkkkyhxhyxyxy 21! 2)( 211() , ()kkkkkkkhy xy xhy xy xy xhfxy x 略略去去 1, 1 , 0,10nkyxhfyyyaykkkk用用 的近似值的近似值 代入上式右端，记所得结果代入上式右端，记所得结果为为，则得到数值解序列的计算公式，则得到数值解序列的计算公式：()ky xky1ky 2、化导数为差商的求解方法思路：、化导数为差商的求解方法思路：若在点若在点 处的导数用差商来近似代替，如向前差商处的导数用差商来近似代替，如向前差商 hxyxyxykkk

6、1kx则微分方程初值问题化为则微分方程初值问题化为 1, 1 , 001nkyayxyhxyxykkk将近似号改为等号，精确解将近似号改为等号，精确解 改为近似解改为近似解 ，得，得ky kxy 1, 1 , 0,10nkyxhfyyyaykkkk3、数值积分的求解思路、数值积分的求解思路：如果将微分方程如果将微分方程 在各小区间在各小区间 上对其两边进行积分，即上对其两边进行积分，即 yxfy, 1, kkxx 111, 1 , 0,kkkkxxxxnkdxxyxfdxy 011,yaydxxyxfxyxykkxxkk如用矩形数值积分公式可得：如用矩形数值积分公式可得： 1, 1 , 0,1

7、0nkyxhfyyyaykkkk以上三种方法推导出同一个数值求解公式以上三种方法推导出同一个数值求解公式: :这个数值公式称为这个数值公式称为欧拉欧拉(Euler)(Euler)公式。公式。 1, 1 , 0,10nkyxhfyyyaykkkk6.1 欧拉方法欧拉方法一、一、 欧拉格式：欧拉格式：x0 x1向前差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为欧拉公式几何意义欧拉公式几何意义 用一条通过初始点的折用一条通过初始点的折线近似表示解曲线线近似表示解曲线 ,亦称为亦称为欧欧拉折线法拉折线法 ，或称为或称为

8、矩形法。矩形法。)1,., 0(),(1 nkyxfhyykkkk一般形式一般形式1 1、显式欧拉公式、显式欧拉公式在假设在假设 yk = y(xk)，即第，即第 k 步计算是精确的前提下，考步计算是精确的前提下，考虑的截断误差虑的截断误差 Rk = y(xk+1) yk+1 称为称为局部截断误差局部截断误差 。定义定义 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。定义定义 欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()(2112kkyxhfyyxyhxyyxyRkhkkkki )(22 yh 欧拉法具有欧拉法

9、具有 1 阶精度。阶精度。局部截断误差和阶数局部截断误差和阶数1 kkxx 2、隐式欧拉格式、隐式欧拉格式向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 nkyxfhyykiik由于未知数由于未知数 yk+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：

10、11)( kkkyxyR)(22 yh 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。二、两步欧拉格式（中点公式）二、两步欧拉格式（中点公式）中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 nkyxfhyykkkk假设假设 ,则可以导出则可以导出即两步欧拉格式具有即两步欧拉格式具有 2 阶精度。阶精度。)(),(11kkkkxyyxyy 311()()kkkRy xyO h 该方法需要该方法需要2个初值个初值 y0和和 y1来启动递推过程，这样的算法称来启动递推过程，这样的算法

11、称为为双步法双步法。三、三、 梯形公式梯形公式 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 nkyxfyxfhyykkkkkk注：注：有局部截断误差有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法迭代法，不易求解。，不易求解。)()(311hOyxyRkkk 对欧拉法进行改进，用梯形公式计算右侧积分，即对欧拉法进行改进，用梯形公式计算右侧积分，即 11,2,1 kkkkxxyxfyxfhdxyxfkk 2 ,

12、 1 , 0,),(),(2),()0(11)1(1)0(1 ikyxfyxfhyyyxhfyykkkkkikkkkk计算计算公式公式梯形格式算法计算步骤：梯形格式算法计算步骤：先用先用(1)式计算出式计算出 处处 。1kx (0)1ky 再用再用(2)式反复进行迭代，得到式反复进行迭代，得到(1)(2)11,kkyy(0)1(1)(0)111(,)(,)(,)2kkkkikkkkkkyyhf xyhyyf xyf xy 计算计算公式公式-(1)-(2)类似地得到类似地得到 用用 控制迭代次数，控制迭代次数， 为允许误差。为允许误差。把满足误差要求的把满足误差要求的 作为作为 的近似值。的近似

13、值。 (1)( )11iikkyy (1)1iky 1ky x 23,kkyy方方 法法 显式欧拉法显式欧拉法 隐式欧拉法隐式欧拉法 梯形公式梯形公式 中点公式中点公式 简单简单 精度低精度低 稳定性最好稳定性最好 精度低精度低, 计算量大计算量大 精度提高精度提高 计算量大计算量大 精度提高精度提高, 显式显式 多一个初值多一个初值, 可能影响精度可能影响精度 不同方法比较不同方法比较四、四、改进欧拉法（预报改进欧拉法（预报-校正法）校正法）Step 1: 先用先用显式显式欧拉公式作欧拉公式作预报预报，算出，算出),(1iiiiyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式

14、的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii它可表示为嵌套形式它可表示为嵌套形式表示为平均化形式表示为平均化形式 pcipiiciiipyyyyxhfyyyxhfyy21,11此法称为此法称为预报预报-校正法，校正法，是是显式算法显式算法。注：注：可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式（只迭代一次）递推格式（只迭代一次） ，比隐式梯形公式的迭代，比隐式梯形公式的迭代求解过程求解过

15、程简单简单。它的。它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。脚标用脚标用 i举例：进行比较。进行比较。精确解精确解并与并与法求解法求解法和改进法和改进试分别用试分别用设初值问题设初值问题例例xyyyxydxdy21,EulerEuler1)0(2 计算结果如下表所示：计算结果如下表所示：法：法：改进的改进的法：法：上结果，此时上结果，此时计算计算解：取解：取 ,.)2 , 1 , 0()(21)2()2(Euler,.)2 , 1 , 0()2( 1 . 01 , 0, 1 . 0111iyyyyxyhyyyxyhyyiyxyyyEulerxhcpipipiciiiipiiiiixEul

16、er法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.73

17、78671.7320516.2 龙格龙格 - 库塔法库塔法一、一、 泰勒级数法泰勒级数法 龙格龙格库塔库塔(Runge-Kut ta)法法(简称为简称为R-K方法方法)是一类高精是一类高精度的一步法，这类方法与泰勒级数法有着密度的一步法，这类方法与泰勒级数法有着密 切的关系。切的关系。 设有初值问题设有初值问题 00)()(,()(yxyxyxfxy由由 泰勒展开式泰勒展开式 1)(21! 2)( kkkkkkkkhxykhxyhxhyxyxy 从理论上讲，只要解从理论上讲，只要解y(x)有任意阶导数，泰勒展开方有任意阶导数，泰勒展开方法就可以构造法就可以构造任意阶任意阶求求yk+1公式。但由

18、于计算这些导数是公式。但由于计算这些导数是非常复杂的，所以这种方法实际上不能用来解初值问题。非常复杂的，所以这种方法实际上不能用来解初值问题。 设有初值问题设有初值问题 00)()(,()(yxyxyxfxy二、二、龙格库塔法的基本思路龙格库塔法的基本思路)(),()(hxyhxhfxykkk 1)(,()()(1kkxxkkdxxyxfxyxy等价于：等价于：（积分中值定理）（积分中值定理） R-K方法基本思想：方法基本思想：用用 在几个不同点的加权平均值在几个不同点的加权平均值（线性组合）来代替准确的（线性组合）来代替准确的 的值，构的值，构造近似公式。再把近似公式与解的泰勒展开造近似公式

19、。再把近似公式与解的泰勒展开 式进行比较，式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。使前面的若干项相同，从而使近似公式达到一定的阶数。 这样龙格库塔法保留了泰勒级数展开法的高阶局部截这样龙格库塔法保留了泰勒级数展开法的高阶局部截断误差，又避免了高阶导数的计算。断误差，又避免了高阶导数的计算。 我们先分析欧拉法我们先分析欧拉法 与预估与预估校正法。校正法。),(yxf)(,(hxyhxfkk 112111121(,)(,)(,)(,)1()()22pkkkkkckkpkkkkkpcyyhf xyKf xyyyhf xyKf xyhKhyyKKyyy 改改进进的的欧欧拉拉式式在

20、在公公中中1111131(,)(,)(,)(,),()kkkkkkkkkkkkxyyf xyxyhKyf xyhKxxO h 则则是是用用点点处处的的斜斜率率和和由由此此点点处处信信息息预预估估的的点点处处的的斜斜率率的的算算术术平平均均值值来来近近似似代代替替区区间间上上的的平平均均斜斜率率，局局部部截截断断误误差差为为。0121( )(,)(,)(,),()kkkkkkkkkkky ayyyhf xyxyyf xyxxO h 在在欧欧拉拉公公式式中中仅仅用用一一个个点点处处的的斜斜率率来来近近似似代代替替区区间间上上的的平平均均斜斜率率，局局部部截截断断误误差差为为。所所以以如如果果在在区

21、区间间上上多多预预估估几几个个点点的的斜斜率率值值，再再将将它它们们的的线线性性组组合合作作为为平平均均斜斜率率的的近近似似值值，则则就就有有可可能能构构造造出出精精度度更更高高的的计计算算格格式式。推广推广1(, )kkkkyyhxyh ( , )f x y 其其中中 是是用用在在一一些些点点上上值值的的线线性性组组合合来来构构成成这种单步法称为这种单步法称为Runge-KuttaRunge-Kutta方法方法, ,简记为简记为R-KR-K公式公式. .1(, )rkkiiixyhc k K Ki i为某些点上的斜率，或为某些点上的斜率，或f(x,y)f(x,y)在某些点上的值。在某些点上的

22、值。,Runge-Kutta.若若 是是由由 个个 值值线线性性组组合合构构成成 则则称称线线性性阶阶方方法法RRf 三、三、 二阶龙格二阶龙格 - 库塔法库塔法目标目标：建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyx

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

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

25、iiiiiii 要求要求 ，则必须有：，则必须有：)()(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 ) 均为待定系数，确定

26、这些系数的步骤与前面相似。其均为待定系数，确定这些系数的步骤与前面相似。其解不唯一。解不唯一。).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四阶最常用为四阶经典龙格经典龙格-库塔法库塔法四阶经典龙格四阶经典龙格-库塔法公式库塔法公式),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 四、四、 四阶龙格四阶龙格 - 库塔法库塔法用四个用四个f

27、 f函数值的线性组合得到四阶函数值的线性组合得到四阶龙格龙格 - - 库塔法库塔法。经典龙格经典龙格-库塔法公式具有四阶精度，因此可取大步长。库塔法公式具有四阶精度，因此可取大步长。注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数2345642n 8n高于四阶时每步计算量增加较多，但精度提高不快，高于四阶时每步计算量增加较多，但精度提

28、高不快，因此使用的比较少。因此使用的比较少。由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。写出经典龙格写出经典龙格- -库塔法求解初值问题库塔法求解初值问题 1)0(1yxyy 的计算公式，并取步长的计算公式，并取步长2 . 0h，计算，计算)4 . 0(y的近似值。的近似值。 解解 由题意知由题意知1,xyyxf，所以有，所以有例 121132211234,1(,)() () 122220.20.21

29、10.90.9122(,)() () 122220.20.20.90.911022226nnnnnnnnnnnnnnnnnnnnnnnnkf x yyxhhhhkf xykykxyyxxyxhhhhkf xykykxyyxxhyykkkk 433.910.911,10.2 0.910.9110.2 10.8180.8181nnnnnnnnnnnnyxkf xh yhkyhkxhyyxxyx 1)0(1yxyy代入，有代入，有 11234220.8187330.1812670.26nnnnhyykkkkyx 由于由于1)0(0 yy，取，取2 . 0, 0010hxxx，取步长，取步长2 . 0

30、h， 有有 018733. 1)2 . 0(1 yy, , 070324. 1)4 . 0(2 yy Lab 15. Runge-Kutta (Order Four) Use Runge-Kutta method of order four to approximate the solution of the initial-value problem, , and (1)You are supposed to write a functionvoid RK4 ( double (*f)( ), double a, double b, double y0, int n, FILE *outfi

31、le)to approximate the solution of Problem (1) with y= f, x in a, b, and the initial value of y being y0. Output the approximating values of y on the n+1 equally spaced grid points from a to b to outfile.InputThere is no input file. Instead, you must hand in your function in a *.h file. The rule of n

32、aming the *.h file is the same as that of naming the *.c or *.cpp files.Output ( represents a space) For each test case, you are supposed to print n+1 lines, and each line is in the following format: fprintf( outfile, “%8.4f%12.8en”, x, y );),(yxfy ,bax 0)(yay Sample Judge ProgramSample Judge Progra

33、m#include #include #include 98115001_15.h double f0(double x, double y) return (y x*x+1.0); void main( ) FILE *outfile = fopen(out.txt, w);int n;double a, 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 ( ( repr

34、esents a space) represents a space) 0.00005.00000000e 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 a

35、nd Stability */ 收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 0)1 (yhyii 对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihy

36、hyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656