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

《常微分方程的数值解法.ppt》由会员分享，可在线阅读，更多相关《常微分方程的数值解法.ppt（33页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程的数值解法1现在学习的是第1页，共33页上一页上一页 下一页下一页 返回返回 本章介绍求解微分方程数值解的基本思想和方法本章介绍求解微分方程数值解的基本思想和方法.含有自变量、未知函数和它的一阶导数和高阶导数的方含有自变量、未知函数和它的一阶导数和高阶导数的方程程.常微分方程常微分方程 它是描述运动、变化规律的重要数学方法之一，分为它是描述运动、变化规律的重要数学方法之一，分为两类：两类：1.初值问题初值问题 即给出未知函数及导数在初始点的值；即给出未知函数及导数在初始点的值；2.边值问题边值问题 即给出未知函数及（或）它的某些导数在区即给出未知函数及（或）它的某些导数在区间两个端点

2、的值间两个端点的值。2现在学习的是第2页，共33页 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题:只要只要 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.节

3、点间距节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi=h(常数常数)。yn 称为问题的称为问题的数值解数值解.数值解所满足的离散方程统称为数值解所满足的离散方程统称为差分格式差分格式.上一页上一页 下一页下一页 返回返回 3现在学习的是第3页，共33页第一节第一节 欧拉方法欧拉方法一、一、欧拉公式欧拉公式令令yn为为y(xn)的近似值，将上式代入的近似值，将上式代入(*)式可得式可得此式称为此式称为欧拉欧拉(Euler)公式公式.为为Euler方法的方法的局部截断误差局部截断误差.上一页上一页 下一页下一页 返回返回 4现在学习的是第4页，共33页例例1 用欧拉公

4、式解初值问题用欧拉公式解初值问题解解：取步长取步长 h=0.1,欧拉公式的具体形式为：欧拉公式的具体形式为：依次计算可得依次计算可得 y3，y4，y5，y6，y7，y8，y9，y10上一页上一页 下一页下一页 返回返回 5现在学习的是第5页，共33页其部分结果见下表其部分结果见下表 可见可见Euler方法的计算结果精度不太高。方法的计算结果精度不太高。上一页上一页 下一页下一页 返回返回 6现在学习的是第6页，共33页 欧拉公式的几何意义：欧拉公式的几何意义：x0P0 x1P1x2P2xnPn几何意义：几何意义：用折线近似代替方程的解曲线用折线近似代替方程的解曲线,因而也称因而也称Euler方

5、法为方法为折线法折线法.上一页上一页 下一页下一页 返回返回 7现在学习的是第7页，共33页二、二、后退的后退的欧拉公式欧拉公式也用一阶差商逼近导数也用一阶差商逼近导数令令yn+1为为y(xn+1)的近似值，则可得的近似值，则可得称为称为后退后退Euler公式公式 已知已知 yn时时,必须通过解方程才能求出必须通过解方程才能求出yn1,这样的公式称为这样的公式称为隐式公隐式公式式,而而Euler公式为公式为显式公式显式公式.Euler公式和后退公式和后退Euler公式都是由公式都是由yn去计算去计算yn+1，因此，称，因此，称它们为它们为单步法。单步法。上一页上一页 下一页下一页 返回返回 8

6、现在学习的是第8页，共33页在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的截断误差的截断误差 Ti+1=y(xi+1)yi+1 称为称为局部截断误差。局部截断误差。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。显然显然,p越大越大,精度越高精度越高.三、三、局部截断误差与方法的阶局部截断误差与方法的阶（将准确解（将准确解代入公式的左、右两端，其左端与右端之差）代入公式的左、右两端，其左端与右端之差）Euler方法的精度方法的精度 其中：其中：上一页上一页 下一页下一页 返回返回 9

7、现在学习的是第9页，共33页所以，所以，Euler方法具有方法具有 1 阶精度。阶精度。将将在点在点处一阶处一阶Taylor展开展开上一页上一页 下一页下一页 返回返回 10现在学习的是第10页，共33页所以，所以，后退的后退的Euler方法也具有方法也具有 1 阶精度。阶精度。将将在点在点处一阶处一阶Taylor展开展开 隐式隐式Euler方法的精度方法的精度 上一页上一页 下一页下一页 返回返回 11现在学习的是第11页，共33页 显、隐式两种算法的显、隐式两种算法的平均平均 欧拉公式的改进欧拉公式的改进其局部误差为：其局部误差为：此公式具有此公式具有2 2阶精度阶精度.称称平均公式或梯形

8、公式平均公式或梯形公式梯形公式可由下迭代式计算：其中迭代初值是梯形公式可由下迭代式计算：其中迭代初值是Euler公式提供公式提供.上一页上一页 下一页下一页 返回返回 12现在学习的是第12页，共33页四、四、改进的欧拉公式改进的欧拉公式Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法。可以证明该算法具有可以证明该算法具有 2 阶精度，同

9、时可以看到它是个阶精度，同时可以看到它是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。它的。它的精度高精度高于显式欧拉法。于显式欧拉法。上一页上一页 下一页下一页 返回返回 13现在学习的是第13页，共33页为了便于编程为了便于编程,常将常将改进的欧拉公式改进的欧拉公式写为：写为：上一页上一页 下一页下一页 返回返回 14现在学习的是第14页，共33页例例2 用用改进的欧拉法解例改进的欧拉法解例1中的初值问题中的初值问题.解：解：取步长取步长 h=0.1,改进欧拉法的具体改进欧拉法的具体 形式为形式为具体计算过程如下具体计算过程如下上一页上一页 下一页

10、下一页 返回返回 15现在学习的是第15页，共33页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现在学习的是第16页，共33页例例3 对下面的对下面的初值问题初值问题解解 （1）取步长取步长 h=0.1,欧拉方法的具体公式为欧

11、拉方法的具体公式为（2）取步长取步长 h=0.1,改进的欧拉方法的具体公式为改进的欧拉方法的具体公式为取步长取步长h=0.1，分别用，分别用Euler方法、改进的方法、改进的Euler方法求数值解。方法求数值解。上一页上一页 下一页下一页 返回返回 17现在学习的是第17页，共33页计算结果见下表计算结果见下表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.3

12、48 6790.905 0000.819 0250.741 2180.670 8020.607 0760.549 4040.497 2100.449 9750.407 2280.368 541上一页上一页 下一页下一页 返回返回 18现在学习的是第18页，共33页第二节第二节 龙格龙格-库塔法库塔法基本思想基本思想考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？只要能对平均斜率提供一种近似算法只要能对平均斜率提供一种近似算法,就能得到一种对应的差分格式就能得到一种对应的差分格式.

13、上一页上一页 下一页下一页 返回返回 19现在学习的是第19页，共33页例如取例如取 m 个点的斜率构造如下形式的公式个点的斜率构造如下形式的公式该公式称为该公式称为m级级龙格库塔龙格库塔(Runge-Kutta)公式公式,简称简称R-K公式公式.求解：求解：只需将公式的局部截断误差在只需将公式的局部截断误差在xn点进行点进行Taylor展开展开,令其前面令其前面尽可能多的项为尽可能多的项为0,便可导出便可导出ai，bij，ci所满足的方程组所满足的方程组,即可从中求出这些即可从中求出这些系数系数.上一页上一页 下一页下一页 返回返回 20现在学习的是第20页，共33页以以 m=2 的情形为例

14、说明建立的情形为例说明建立R-K公式的方法公式的方法.其局部截断误差为：其局部截断误差为：上一页上一页 下一页下一页 返回返回 21现在学习的是第21页，共33页因此有：因此有：而对于而对于h3，若将若将k2的的Taylor展开式多取一项展开式多取一项,会发现会发现h3项的系数不可项的系数不可能为能为0.而对于上式有无穷多个解而对于上式有无穷多个解,它的每一组解都给出了一个局部截它的每一组解都给出了一个局部截断误差为断误差为 的二级的二级R-K公式公式,即即二阶二阶R-K公式公式.当取当取时时,二阶二阶R-K公式就是改进的公式就是改进的Euler公式公式 这里有这里有 个未知数，个未知数，个方

15、程。个方程。32上一页上一页 下一页下一页 返回返回 22现在学习的是第22页，共33页 常用的标准四阶常用的标准四阶RK公式公式（经典（经典R-K方法）方法）最常用的四阶标准最常用的四阶标准RK公式（经典公式（经典R-K方法）为：方法）为：上一页上一页 下一页下一页 返回返回 23现在学习的是第23页，共33页例例 用四阶标准用四阶标准R-K公式解初值问题公式解初值问题 解：取解：取 h=0.2，四阶标准，四阶标准R-K法的具体格式如下法的具体格式如下:上一页上一页 下一页下一页 返回返回 24现在学习的是第24页，共33页已知已知 上一页上一页 下一页下一页 返回返回 25现在学习的是第2

16、5页，共33页同理可计算得同理可计算得具体结果见下表具体结果见下表 至少具有四位有效数字至少具有四位有效数字.比较比较:上节用改进的上节用改进的Euler公式计算公式计算,取取h=0.1，最多具有四位有效数字，最多具有四位有效数字。上一页上一页 下一页下一页 返回返回 26现在学习的是第26页，共33页 改进的改进的Euler公式每前进一步只要计算两次公式每前进一步只要计算两次f 值值,而而4 4阶阶R-K公式每前公式每前进一步要计算四次进一步要计算四次f 值值,但改进的但改进的Euler法的步长比法的步长比4阶阶R-K法的小一半法的小一半,两者计算总量差不多两者计算总量差不多.而而4阶阶R-

17、K法的效果要比改进的法的效果要比改进的Euler法好法好.由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受解库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。上一页上一页 下一页下一页 返回返回 27现在学习的是第27页，共33页第三节第三节 单步法的单步法的收敛性与稳定性收敛性与稳定性 收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定的的 x=xi=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则

18、称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 上一页上一页 下一页下一页 返回返回 28现在学习的是第28页，共33页 稳定性稳定性/*Stability*/例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法隐式欧拉

19、法隐式欧拉法欧拉法欧拉法 节点节点 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现在学习的是第29页，共33页若若某某算算法

20、法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的.一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程常数常数l l00，可，可以是复数以是复数 当步长取为当步长取为 h 时，将某算法应用于上式，并假设在初时，将某算法应用于上式，并假设在初值产生误差值产生误差 ，则若此误差以后逐步衰减，就称该算，则若此误差以后逐步衰减，就称该算法相对于法相对于z=l l h 绝对稳定绝对稳定，z 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定

21、稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。上一页上一页 下一页下一页 返回返回 30现在学习的是第30页，共33页例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。上一页上一页 下一页下一页 返回返回 31现在学习的是第31页，共33页一些单步法的绝对稳定

22、区间见下表一些单步法的绝对稳定区间见下表方方 法法绝对稳定区间绝对稳定区间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现在学习的是第32页，共33页例例 分别取分别取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现在学习的是第33页，共33页