常微分方程的数值解法讲稿.ppt

常微分方程的数值解法1第一页，讲稿共三十三页哦上一页上一页 下一页下一页 返回返回 本章介绍求解微分方程数值解的基本思想和方法本章介绍求解微分方程数值解的基本思想和方法.含有自变量、未知函数和它的一阶导数和高阶导数含有自变量、未知函数和它的一阶导数和高阶导数的方程的方程.常微分方程常微分方程 它是描述运动、变化规律的重要数学方法之一，分为它是描述运动、变化规律的重要数学方法之一，分为两类：两类：1.初值问题初值问题 即给出未知函数及导数在初始点的值；即给出未知函数及导数在初始点的值；2.边值问题边值问题 即给出未知函数及（或）它的某些导数在区间即给出未知函数及（或）它的某些导数在区间两个端点的值两个端点的值。2第二页，讲稿共三十三页哦 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题:只要只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与即存在与 x,y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a,b 上的上的 y1(x)和和 y2(x)都成立，则上述问题解都成立，则上述问题解存在唯一解存在唯一解。所谓所谓数值解法数值解法就是要计算出初值问题的解函数就是要计算出初值问题的解函数 y(x)在一系列在一系列离散点离散点 a=x0 x1 xN=b上的近似值：上的近似值：y0,y1,yN.节点间距节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi=h(常数常数)。yn 称为问题的称为问题的数值解数值解.数值解所满足的离散方程统称为数值解所满足的离散方程统称为差分格式差分格式.上一页上一页 下一页下一页 返回返回 3第三页，讲稿共三十三页哦第一节第一节 欧拉方法欧拉方法一、一、欧拉公式欧拉公式令令yn为为y(xn)的近似值，将上式代入的近似值，将上式代入(*)式可得式可得此式称为此式称为欧拉欧拉(Euler)公式公式.为为Euler方法的方法的局部截断误差局部截断误差.上一页上一页 下一页下一页 返回返回 4第四页，讲稿共三十三页哦例例1 用欧拉公式解初值问题用欧拉公式解初值问题解解：取步长取步长 h=0.1,欧拉公式的具体形式为：欧拉公式的具体形式为：依次计算可得依次计算可得 y3，y4，y5，y6，y7，y8，y9，y10上一页上一页 下一页下一页 返回返回 5第五页，讲稿共三十三页哦其部分结果见下表其部分结果见下表 可见可见Euler方法的计算结果精度不太高。方法的计算结果精度不太高。上一页上一页 下一页下一页 返回返回 6第六页，讲稿共三十三页哦 欧拉公式的几何意义：欧拉公式的几何意义：x0P0 x1P1x2P2xnPn几何意义：几何意义：用折线近似代替方程的解曲线用折线近似代替方程的解曲线,因而也称因而也称Euler方法为方法为折线法折线法.上一页上一页 下一页下一页 返回返回 7第七页，讲稿共三十三页哦二、二、后退的后退的欧拉公式欧拉公式也用一阶差商逼近导数也用一阶差商逼近导数令令yn+1为为y(xn+1)的近似值，则可得的近似值，则可得称为称为后退后退Euler公式公式 已知已知 yn时时,必须通过解方程才能求出必须通过解方程才能求出yn1,这样的公式称为这样的公式称为隐式公式隐式公式,而而Euler公式为公式为显式公式显式公式.Euler公式和后退公式和后退Euler公式都是由公式都是由yn去计算去计算yn+1，因此，称，因此，称它们为它们为单步法。单步法。上一页上一页 下一页下一页 返回返回 8第八页，讲稿共三十三页哦在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的截断误差的截断误差 Ti+1=y(xi+1)yi+1 称为称为局部截断误差。局部截断误差。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。显然显然,p越大越大,精度越高精度越高.三、三、局部截断误差与方法的阶局部截断误差与方法的阶（将准确解（将准确解代入公式的左、右两端，其左端与右端之差）代入公式的左、右两端，其左端与右端之差）Euler方法的精度方法的精度 其中：其中：上一页上一页 下一页下一页 返回返回 9第九页，讲稿共三十三页哦所以，所以，Euler方法具有方法具有 1 阶精度。阶精度。将将在点在点处一阶处一阶Taylor展开展开上一页上一页 下一页下一页 返回返回 10第十页，讲稿共三十三页哦所以，所以，后退的后退的Euler方法也具有方法也具有 1 阶精度。阶精度。将将在点在点处一阶处一阶Taylor展开展开 隐式隐式Euler方法的精度方法的精度 上一页上一页 下一页下一页 返回返回 11第十一页，讲稿共三十三页哦 显、隐式两种算法的显、隐式两种算法的平均平均 欧拉公式的改进欧拉公式的改进其局部误差为：其局部误差为：此公式具有此公式具有2 2阶精度阶精度.称称平均公式或梯形公式平均公式或梯形公式梯形公式可由下迭代式计算：其中迭代初值是梯形公式可由下迭代式计算：其中迭代初值是Euler公式提供公式提供.上一页上一页 下一页下一页 返回返回 12第十二页，讲稿共三十三页哦四、四、改进的欧拉公式改进的欧拉公式Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。它的。它的精精度高度高于显式欧拉法。于显式欧拉法。上一页上一页 下一页下一页 返回返回 13第十三页，讲稿共三十三页哦为了便于编程为了便于编程,常将常将改进的欧拉公式改进的欧拉公式写为：写为：上一页上一页 下一页下一页 返回返回 14第十四页，讲稿共三十三页哦例例2 用用改进的欧拉法解例改进的欧拉法解例1中的初值问题中的初值问题.解：解：取步长取步长 h=0.1,改进欧拉法的具体改进欧拉法的具体 形式为形式为具体计算过程如下具体计算过程如下上一页上一页 下一页下一页 返回返回 15第十五页，讲稿共三十三页哦xn改进的欧拉法改进的欧拉法误差误差xn改进的欧拉法改进的欧拉法误差误差0100.61.4859560.0027160.21.1840960.0000880.81.6164760.0040240.41.3433600.0017191.01.7378690.005818依次计算可得依次计算可得 y3，y4，y5，y6，y7，y8，y9，y10其部分结果见下表其部分结果见下表 上一页上一页 下一页下一页 返回返回 16第十六页，讲稿共三十三页哦例例3 对下面的对下面的初值问题初值问题解解 （1）取步长取步长 h=0.1,欧拉方法的具体公式为欧拉方法的具体公式为（2）取步长取步长 h=0.1,改进的欧拉方法的具体公式为改进的欧拉方法的具体公式为取步长取步长h=0.1，分别用，分别用Euler方法、改进的方法、改进的Euler方法求数值解。方法求数值解。上一页上一页 下一页下一页 返回返回 17第十七页，讲稿共三十三页哦计算结果见下表计算结果见下表Euler方法方法改进的改进的Euler方法方法xnynyn0.10.20.30.40.50.60.70.80.91.00.900 0000.810 0000.729 0000.656 1000.590 4900.531 4410.478 2970.430 4670.387 4210.348 6790.905 0000.819 0250.741 2180.670 8020.607 0760.549 4040.497 2100.449 9750.407 2280.368 541上一页上一页 下一页下一页 返回返回 18第十八页，讲稿共三十三页哦第二节第二节 龙格龙格-库塔法库塔法基本思想基本思想考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？只要能对平均斜率提供一种近似算法只要能对平均斜率提供一种近似算法,就能得到一种对应的差分格式就能得到一种对应的差分格式.上一页上一页 下一页下一页 返回返回 19第十九页，讲稿共三十三页哦例如取例如取 m 个点的斜率构造如下形式的公式个点的斜率构造如下形式的公式该公式称为该公式称为m级级龙格库塔龙格库塔(Runge-Kutta)公式公式,简称简称R-K公式公式.求解：求解：只需将公式的局部截断误差在只需将公式的局部截断误差在xn点进行点进行Taylor展开展开,令其前面尽令其前面尽可能多的项为可能多的项为0,便可导出便可导出ai，bij，ci所满足的方程组所满足的方程组,即可从中求出这些即可从中求出这些系数系数.上一页上一页 下一页下一页 返回返回 20第二十页，讲稿共三十三页哦以以 m=2 的情形为例说明建立的情形为例说明建立R-K公式的方法公式的方法.其局部截断误差为：其局部截断误差为：上一页上一页 下一页下一页 返回返回 21第二十一页，讲稿共三十三页哦因此有：因此有：而对于而对于h3，若将若将k2的的Taylor展开式多取一项展开式多取一项,会发现会发现h3项的系数不可能项的系数不可能为为0.而对于上式有无穷多个解而对于上式有无穷多个解,它的每一组解都给出了一个局部截它的每一组解都给出了一个局部截断误差为断误差为 的二级的二级R-K公式公式,即即二阶二阶R-K公式公式.当取当取时时,二阶二阶R-K公式就是改进的公式就是改进的Euler公式公式 这里有这里有 个未知数，个未知数，个方程。个方程。32上一页上一页 下一页下一页 返回返回 22第二十二页，讲稿共三十三页哦 常用的标准四阶常用的标准四阶RK公式公式（经典（经典R-K方法）方法）最常用的四阶标准最常用的四阶标准RK公式（经典公式（经典R-K方法）为：方法）为：上一页上一页 下一页下一页 返回返回 23第二十三页，讲稿共三十三页哦例例 用四阶标准用四阶标准R-K公式解初值问题公式解初值问题 解：取解：取 h=0.2，四阶标准，四阶标准R-K法的具体格式如下法的具体格式如下:上一页上一页 下一页下一页 返回返回 24第二十四页，讲稿共三十三页哦已知已知 上一页上一页 下一页下一页 返回返回 25第二十五页，讲稿共三十三页哦同理可计算得同理可计算得具体结果见下表具体结果见下表 至少具有四位有效数字至少具有四位有效数字.比较比较:上节用改进的上节用改进的Euler公式计算公式计算,取取h=0.1，最多具有四位有效数字，最多具有四位有效数字。上一页上一页 下一页下一页 返回返回 26第二十六页，讲稿共三十三页哦 改进的改进的Euler公式每前进一步只要计算两次公式每前进一步只要计算两次f 值值,而而4 4阶阶R-K公式公式每前进一步要计算四次每前进一步要计算四次f 值值,但改进的但改进的Euler法的步长比法的步长比4阶阶R-K法的法的小一半小一半,两者计算总量差不多两者计算总量差不多.而而4阶阶R-K法的效果要比改进的法的效果要比改进的Euler法好法好.由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用的光滑性影响。对于光滑性不太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。上一页上一页 下一页下一页 返回返回 27第二十七页，讲稿共三十三页哦第三节第三节 单步法的单步法的收敛性与稳定性收敛性与稳定性 收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定的的 x=xi=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 上一页上一页 下一页下一页 返回返回 28第二十八页，讲稿共三十三页哦 稳定性稳定性/*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 7What is wrong?!上一页上一页 下一页下一页 返回返回 29第二十九页，讲稿共三十三页哦若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的.一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程常数常数l l00，可以是复数可以是复数 当步长取为当步长取为 h 时，将某算法应用于上式，并假设在初值时，将某算法应用于上式，并假设在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该算法，则若此误差以后逐步衰减，就称该算法相对于相对于z=l l h 绝对稳定绝对稳定，z 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。上一页上一页 下一页下一页 返回返回 30第三十页，讲稿共三十三页哦例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。上一页上一页 下一页下一页 返回返回 31第三十一页，讲稿共三十三页哦一些单步法的绝对稳定区间见下表一些单步法的绝对稳定区间见下表方方 法法绝对稳定区间绝对稳定区间Euler方法方法改进的改进的Euler方法方法三阶三阶R-K法法四阶四阶R-K法法隐式隐式Euler法法梯形法梯形法-2 z 0-2 z 0-2.51 z 0-2.785 z 0-z 0-l l 0上一页上一页 下一页下一页 返回返回 32第三十二页，讲稿共三十三页哦例例 分别取分别取h=1，2，4，用经典四阶，用经典四阶R-K方法求解方法求解讨论稳定性。讨论稳定性。解：解：由上表知道，四阶由上表知道，四阶R-K法的绝对稳定区间是法的绝对稳定区间是-2.785 z 0，即：当即：当h 2.785时，四阶时，四阶R-K方法才稳定。方法才稳定。故，当故，当h=1，2时，四阶时，四阶R-K方法的计算结果稳定，方法的计算结果稳定，h=4的时候，的时候，结果不稳定。结果不稳定。上一页上一页 下一页下一页 返回返回 33第三十三页，讲稿共三十三页哦