常微分方程初值问题的数值解法 (2)讲稿.ppt

常微分方程初值问题的数值解法(2)第一页，讲稿共三十五页哦引言引言在实际问题中，常需要求解微分方程在实际问题中，常需要求解微分方程(如发电机转子运动方程如发电机转子运动方程)。只有简单的和典型的微分方程可以求出解析解，而在实只有简单的和典型的微分方程可以求出解析解，而在实际问题中的微分方程往往无法求出解析解。际问题中的微分方程往往无法求出解析解。常微分方程：常微分方程：-(1)-(2)一阶常微分方程一阶常微分方程第二页，讲稿共三十五页哦-(3)(1),(2)式称为式称为初值问题初值问题,(,(3)式称为式称为边值问题边值问题-(4)另外另外,在实际应用中还经常需要求解常微分方程组在实际应用中还经常需要求解常微分方程组:本课程主要研究问题本课程主要研究问题一阶常微分方程一阶常微分方程(1)的数值解法的数值解法,我们首先介绍初值问题我们首先介绍初值问题(1)(1)的解存在的条件的解存在的条件第三页，讲稿共三十五页哦定理定理 只要只要 f(x,y)连续，且关于连续，且关于 y 满足满足 Lipschitz 条件条件，即，即存在与存在与 x,y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a,b 上的上的 y1(x)和和 y2(x)都成立，则初值问题都成立，则初值问题（1）存在唯一解存在唯一解。（通常采用（通常采用等距节点）等距节点）对于问题对于问题(1)(1)要求它的要求它的数值解数值解第四页，讲稿共三十五页哦常微分方程数值解公式的推导常微分方程数值解公式的推导 求初值问题数值解的方法是求初值问题数值解的方法是步进法步进法，即从已知的初值，即从已知的初值y0出出发，通过发，通过一定的计算一定的计算求求y1，然后由，然后由y1或或y0和和y1求出求出y2，依次计，依次计算到算到yn，即在计算出，即在计算出yk后计算后计算yk+1，这时有，这时有单步法单步法：计算：计算yk+1时，只利用时，只利用yk多步法多步法：计算：计算yk+1时，用到时，用到yk,yk-1,yk-2,常微分方程数值解公式的常微分方程数值解公式的主要推导方法主要推导方法泰勒展开泰勒展开利用差商利用差商利用数值积分法利用数值积分法第五页，讲稿共三十五页哦1、泰勒展开的求解思路、泰勒展开的求解思路：将将 按泰勒级数展开按泰勒级数展开用用 的近似值的近似值 代入上式右端，记所得结果代入上式右端，记所得结果为为，则得到数值解序列的计算公式，则得到数值解序列的计算公式：第六页，讲稿共三十五页哦2、化导数为差商的求解方法思路：、化导数为差商的求解方法思路：若在点若在点 处的导数用差商来近似代替，如向前差商处的导数用差商来近似代替，如向前差商则微分方程初值问题化为则微分方程初值问题化为将近似号改为等号，精确解将近似号改为等号，精确解 改为近似解改为近似解 ，得，得第七页，讲稿共三十五页哦3、数值积分的求解思路、数值积分的求解思路：如果将微分方程如果将微分方程 在各小区间在各小区间 上对其两边进行积分，即上对其两边进行积分，即如用矩形数值积分公式可得：如用矩形数值积分公式可得：第八页，讲稿共三十五页哦以上三种方法推导出同一个数值求解公式以上三种方法推导出同一个数值求解公式:这个数值公式称为这个数值公式称为欧拉欧拉(Euler)(Euler)公式。公式。第九页，讲稿共三十五页哦6.1 欧拉方法欧拉方法一、一、欧拉格式：欧拉格式：x0 x1向前差商近似导数向前差商近似导数记为记为欧拉公式几何意义欧拉公式几何意义 用一条通过初始点的折线近用一条通过初始点的折线近似表示解曲线似表示解曲线,亦称为亦称为欧拉折线法欧拉折线法，或称为或称为矩形法。矩形法。一般形式一般形式1 1、显式欧拉公式、显式欧拉公式第十页，讲稿共三十五页哦在在假假设设 yk=y(xk)，即即第第 k 步步计计算算是是精精确确的的前前提提下下，考考虑的截断误差虑的截断误差 Rk=y(xk+1)yk+1 称为称为局部截断误差局部截断误差。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。欧拉法的局部截断误差：欧拉法的局部截断误差：欧拉法具有欧拉法具有 1 阶精度。阶精度。局部截断误差和阶数局部截断误差和阶数第十二页，讲稿共三十五页哦2、隐式欧拉格式、隐式欧拉格式向后差商近似导数向后差商近似导数x0 x1)(,()(1101xyxfhyxy+)1,.,0(),(111=+=+nkyxfhyykiik由于未知数由于未知数 yk+1 同时出现在等式的两边，不能直接得到，故称同时出现在等式的两边，不能直接得到，故称为为隐式隐式 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。第十三页，讲稿共三十五页哦二、两步欧拉格式（中点公式）二、两步欧拉格式（中点公式）中心差商近似导数中心差商近似导数x0 x2x1假设假设 ,则可以导出则可以导出即两步欧拉格式具有即两步欧拉格式具有 2 阶精度。阶精度。该方法需要该方法需要2个初值个初值 y0和和 y1来启动递推过程，这样的算法称为来启动递推过程，这样的算法称为双双步法步法。第十四页，讲稿共三十五页哦三、三、梯形公式梯形公式 显、隐式两种算法的显、隐式两种算法的平均平均注：注：有局部截断误差有局部截断误差 ，即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法迭代法，不易求解。，不易求解。对欧拉法进行改进，用梯形公式计算右侧积分，即对欧拉法进行改进，用梯形公式计算右侧积分，即计算计算公式公式第十五页，讲稿共三十五页哦梯形格式算法计算步骤：梯形格式算法计算步骤：先用先用(1)式计算出式计算出 处处 。再用再用(2)式反复进行迭代，得到式反复进行迭代，得到计算计算公式公式-(1)-(2)类似地得到类似地得到 用用 控制迭代次数，控制迭代次数，为允许误差。把为允许误差。把满足误差要求的满足误差要求的 作为作为 的近似值。的近似值。第十六页，讲稿共三十五页哦方方 法法 显式欧拉法显式欧拉法 隐式欧拉法隐式欧拉法 梯形公式梯形公式 中点公式中点公式 简单简单 精度低精度低 稳定性最好稳定性最好 精度低精度低,计算量大计算量大 精度提高精度提高 计算量大计算量大 精度提高精度提高,显式显式 多一个初值多一个初值,可能影响精度可能影响精度 不同方法比较不同方法比较第十七页，讲稿共三十五页哦四、四、改进欧拉法（预报改进欧拉法（预报-校正法）校正法）Step 1:先用先用显式显式欧拉公式作欧拉公式作预报预报，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy它可表示为嵌套形式它可表示为嵌套形式表示为平均化形式表示为平均化形式此法称为此法称为预报预报-校正法，校正法，是是显式算法显式算法。注：注：可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式（只迭代一次）递推格式（只迭代一次），比隐式梯形公式的迭代求，比隐式梯形公式的迭代求解过程解过程简单简单。脚标用脚标用 i第十八页，讲稿共三十五页哦举例：第十九页，讲稿共三十五页哦第二十页，讲稿共三十五页哦xEuler法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.7378671.732051第二十一页，讲稿共三十五页哦6.2 龙格龙格-库塔法库塔法一、一、泰勒级数法泰勒级数法 龙格龙格库塔库塔(Runge-Kut ta)法法(简称为简称为R-K方法方法)是一类高精度的一步是一类高精度的一步法，这类方法与泰勒级数法有着密法，这类方法与泰勒级数法有着密 切的关系。切的关系。设有初值问题设有初值问题由由 泰勒展开式泰勒展开式 从理论上讲，只要解从理论上讲，只要解y(x)有任意阶导数，泰勒展开方法就可有任意阶导数，泰勒展开方法就可以构造以构造任意阶任意阶求求yk+1公式。但由于计算这些导数是非常复杂的，所公式。但由于计算这些导数是非常复杂的，所以这种方法实际上不能用来解初值问题。以这种方法实际上不能用来解初值问题。第二十二页，讲稿共三十五页哦设有初值问题设有初值问题二、二、龙格库塔法的基本思路龙格库塔法的基本思路等价于：等价于：（积分中值定理）（积分中值定理）R-K方法基本思想：方法基本思想：用用 在几个不同点的加权平均值在几个不同点的加权平均值（线性组合）来代替准确的（线性组合）来代替准确的 的值，构造近的值，构造近似公式。再把近似公式与解的泰勒展开似公式。再把近似公式与解的泰勒展开 式进行比较，使前面式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。的若干项相同，从而使近似公式达到一定的阶数。这样龙格库塔法保留了泰勒级数展开法的高阶局部截这样龙格库塔法保留了泰勒级数展开法的高阶局部截断误差，又避免了高阶导数的计算。断误差，又避免了高阶导数的计算。我们先分析欧拉法我们先分析欧拉法 与预估与预估校正法。校正法。第二十三页，讲稿共三十五页哦第二十四页，讲稿共三十五页哦推广推广这种单步法称为这种单步法称为Runge-KuttaRunge-Kutta方法方法,简记为简记为R-KR-K公式公式.K Ki i为某些点上的斜率，或为某些点上的斜率，或f(x,y)f(x,y)在某些点上的值。在某些点上的值。第二十五页，讲稿共三十五页哦三、三、二阶龙格二阶龙格-库塔法库塔法目标目标：建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能点。欧拉法及其各种变形所能达到的最高精度为达到的最高精度为2阶阶。考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？脚标用脚标用 i第二十六页，讲稿共三十五页哦首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶精阶精度，即在度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入代入yi+1表达式，得到表达式，得到第二十七页，讲稿共三十五页哦Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点的点的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知数，个未知数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库塔格库塔格式式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？第二十八页，讲稿共三十五页哦其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为均为待定系数，确定这些系数的步骤与前面相似。其解不唯一。待定系数，确定这些系数的步骤与前面相似。其解不唯一。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四阶最常用为四阶经典龙格经典龙格-库塔法库塔法第二十九页，讲稿共三十五页哦四阶经典龙格四阶经典龙格-库塔法公式库塔法公式四、四、四阶龙格四阶龙格-库塔法库塔法用四个用四个f f函数值的线性组合得到四阶函数值的线性组合得到四阶龙格龙格 -库塔法库塔法。经典龙格经典龙格-库塔法公式具有四阶精度，因此可取大步长。库塔法公式具有四阶精度，因此可取大步长。第三十页，讲稿共三十五页哦注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的值。的值。Butcher 于于1965年给出了计算量与可达到的最高精度阶年给出了计算量与可达到的最高精度阶数的关系：数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数高于四阶时每步计算量增加较多，但精度提高不快，高于四阶时每步计算量增加较多，但精度提高不快，因此使用的比较少。因此使用的比较少。由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。第三十一页，讲稿共三十五页哦3 收敛性与稳定性收敛性与稳定性 /*Convergency and Stability*/收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定的的 x=xi=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 第三十二页，讲稿共三十五页哦 稳定性稳定性/*Stability*/例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xi 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 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7第三十三页，讲稿共三十五页哦若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的/*absolutely stable*/。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程/*test equation*/常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该算法，则若此误差以后逐步衰减，就称该算法相对于相对于 绝对稳定绝对稳定，的全体构成的全体构成绝对稳定区域绝对稳定区域。我们。我们称称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。h h=h第三十四页，讲稿共三十五页哦例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。第三十五页，讲稿共三十五页哦