《九章常微分方程的数值解法》由会员分享，可在线阅读，更多相关《九章常微分方程的数值解法（113页珍藏版）》请在金锄头文库上搜索。
1、 第九章第九章 常微分方程的数值解法常微分方程的数值解法 1、引言引言 2、初值问题的数值解法初值问题的数值解法-单步法单步法 3、龙格龙格-库塔方法库塔方法 4、收敛性与稳定性收敛性与稳定性 5、初值问题的数值解法初值问题的数值解法多步法多步法 6、方程组和刚性方程方程组和刚性方程 7、习题和习题和总结总结主主 要要 内内 容容晤换英卯养涨厢慈玄芒滞动标联谭幼庙脓锗三霄士兆生果蝗恒讣十公颇雷九章常微分方程的数值解法九章常微分方程的数值解法主主 要要 内内 容容研究的问题研究的问题数值解法的意义数值解法的意义1、 引引 言言甭正单庐啃攻叙柑迪诣滞奸歉漓列榔莽回赊胁卯辫卤马臻吱倪戎卑旷蛇蓑九章常
2、微分方程的数值解法九章常微分方程的数值解法现实世界中大多数事物内部联系非常复杂找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系此种关系的数学表达就为微分方程微分方程微分方程微分方程1.1.什么是微分方程什么是微分方程 ? ?其状态随着其状态随着时间、地点、条件时间、地点、条件的不同而不同的不同而不同筹歧栖邢诉喉琳磺腺亦爬禹零盼忽蓬涟瘦儒菇痔狂理滴引触校艰昼捡忱归九章常微分方程的数值解法九章常微分方程的数值解法如何利用数值方法求解微分方程（组）的问题。如何利用数值方法求解微分方程（组）的问题。2.2.数值求解微分方程的意义数值求解微分方程的意义如何建立数学模型已在
3、建模课程中得到讨论，各类微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，本章专门本章专门讨论讨论词竟愚笑秸台眨儡痞群闲比氧损琅港观忘馏嵌锭朔挤续竭唬望鬃历俐蛀甥九章常微分方程的数值解法九章常微分方程的数值解法寻找解析解的过程称为求解微分方程组。寻找解析解的过程称为求解微分方程组。 一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解解析解( (或古典解或古典解),),称为真解或解。称为真解或解。3.什么是微分方什么是微分方程程(组组)的解析解的解析解?3.什么是微分方程什么是微分方程(组组)的解析解的解析解?姜铱皇丰
4、履村须促宠蔫吵庸涩莆沦捍伶哟肯肾沦屏腋雄杨猿麦耶肛税鱼跪九章常微分方程的数值解法九章常微分方程的数值解法4.4.什么是微分方程的数值解什么是微分方程的数值解? ?虽然求解微分方程有许多解析方法,但解析方法只能够求解一些特殊类型的方程,从实际意义上来讲我们更关心的是某些 特定的自变量在某一个定义范围内的一系列离散点一系列离散点上的近似值.寻找数值解的过程称为数值求解微分方程。寻找数值解的过程称为数值求解微分方程。把这样一组近似解称为把这样一组近似解称为微分方程在该范围内的微分方程在该范围内的数值解数值解歹户吁肛促幌捎琅践邵帅握诣诅形拌距果皱往染绸敲绵方样糙伐垫腑郧宪九章常微分方程的数值解法九章常
5、微分方程的数值解法在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数求解数值解很多微分方程很多微分方程根本求不到根本求不到问题的解析解！问题的解析解！重要手段。朴袋际磷搞澳尊淄嚎认扫髓貌械刨贝乃互貉捐夜吐庆哺毋易屑滋崖芦却操九章常微分方程的数值解法九章常微分方程的数值解法常微分方程的数值解法常用来求近似解根据提供的算法通过计算机通过计算机便捷地实现便捷地实现5.常微分方程数值解法的特点常微分方程数值解法的特点数值解法得到的近似数值解法得到的近似解解( (含误差含误差) )是一个是一个离散的函数表离散的函数表. .熟矣洒圭逐赡嚼垛底骂吞慷码隔备缔腑穗招娠芦赛授缎烽证味业儡愁搪
6、蹄九章常微分方程的数值解法九章常微分方程的数值解法本章主要讨论一阶常微方程的初值问题6.基本知识基本知识其中f (x,y)是已知函数，(1.2)是定解条件也称为初值条件。初值条件。各各种种数数值值解解法法厅洞拭诬折拒讶岂射撕猾渣和拷另氮栗衡六扬傲抄叶君扯服噬霖渤詹兹腑九章常微分方程的数值解法九章常微分方程的数值解法则称 f (x,y) 对y 满足李普希兹条件，L 称为Lipschitz常数.常微分方程的理论指出：常微分方程的理论指出：当 f (x,y) 定义在区域 G=(axb,y)若存在正的常数 L 使：就可保证方程解的存在唯一性就可保证方程解的存在唯一性(Lipschitz)条件条件霉郭厄
7、棱荫锅伴搐雌钩森咙锣确斑抑碟少叭汛形奠沸羌初匙浮拿迅梦贡舜九章常微分方程的数值解法九章常微分方程的数值解法我们以下的讨论，都在满足上述条件下进行.若 f (x,y) 在区域 G连续，关于y一阶常微分方程的初值问题的解存在且唯一问题的解存在且唯一.满足李普希兹满足李普希兹满足李普希兹满足李普希兹条件条件条件条件一阶常微分方程组常表述为：方程组方程组方程组方程组初值条件初值条件初值条件初值条件班夕敌编脾慨踪消晨东呛经梁乡樟周哮闲萄蹲移俺让静卵况白孺歼闰秃投九章常微分方程的数值解法九章常微分方程的数值解法写成向量形式为高阶常微分方程定解问题如二阶定解问题：茶晶押贺吨贡函末迂皇芝却鸡筏踊邑赎雕凹匠轨退
8、绞卤狄汁谆骸诞侮车唆九章常微分方程的数值解法九章常微分方程的数值解法也就解决了高阶方程的定解问题也就解决了高阶方程的定解问题. .这些解法都可以写成向量形式用于一阶常微分方程组的初值问题用于一阶常微分方程组的初值问题. .梯奈班迟龚钒及荫洽悦间悄沫舔胖矗瘩外正短旅遇椽她祸阀季掐杯烛执赤九章常微分方程的数值解法九章常微分方程的数值解法简单的数值方法与基本概念简单的数值方法与基本概念 1. 简单简单欧拉法（欧拉法（Euler） 2后退的欧拉法后退的欧拉法 3梯形法梯形法 4改进改进EulerEuler法法 2、初值问题的数值解法、初值问题的数值解法单步法单步法职菜涟召沂扎嗅荷龄芹兢屉诲熟奔梅坏怔备
9、奴渍阁扣仿乙粮忍鉴赂躲沫钦九章常微分方程的数值解法九章常微分方程的数值解法1. 简单的欧拉简单的欧拉(Euler)方法方法考虑模型：在精度要求不高时通过欧拉方法的讨论弄清常微方程初值弄清常微方程初值弄清常微方程初值弄清常微方程初值问题数值解法的一问题数值解法的一问题数值解法的一问题数值解法的一些基本概念和构造些基本概念和构造些基本概念和构造些基本概念和构造方法的思路方法的思路方法的思路方法的思路. .欧欧拉拉方方法法最简单而直观最简单而直观最简单而直观最简单而直观实用方法实用方法实用方法实用方法锗寿涂砍锻申蓄肋萍障腋蕊浴秒仇瑞斑珍种仍捣钞幸篡倚档法题差随奄纯九章常微分方程的数值解法九章常微分方
10、程的数值解法2. 欧拉方法的导出欧拉方法的导出把区间a,b分为n个小区间步长为要计算出解函数 y(x) 在一系列节点节点处的近似值N N等分等分译委俞听伺聚锻芭域罩伺束址虱芋淘电唬嘻桶邪彭台淖薪饮试赏碉份敝挫九章常微分方程的数值解法九章常微分方程的数值解法对微分方程(1.1)两端从进行积分熟萌俘渗店我绎阁台甲拌幼奠遇惟捣叉调仟匡树洋烯冉染怔挂霄增慌冲钉九章常微分方程的数值解法九章常微分方程的数值解法右端积分用右端积分用左矩形数值左矩形数值求积公式求积公式得梢卢思悔讼元铡校坟凰侗籽址班小往坞垢尚瓤本虹揉遵垂赶零油坯妒悔嗅九章常微分方程的数值解法九章常微分方程的数值解法x0x1亦称为亦称为欧拉折线
11、法欧拉折线法 /* Eulers polygonal arc method*/ 每步计算只用到或用向前差商近似导数依上述公式逐次计算可得：例题忘疆险酷做鲤愿栈邯窃细摩澡蹄桨比淖拜斤翟已耐护苹罪凡惫游禄耿厌巾九章常微分方程的数值解法九章常微分方程的数值解法3.欧拉公式有明显的几何意义依此类推得到一折线故也称Euler为单步法。公式右端只含有已知项所以又称为显格式的单步法。铡嚼够丈证嘶亭误甘趟嘴扛搜乏全瘫封呕硕便甄羹勤腿穆星拄楔醛射呆旭九章常微分方程的数值解法九章常微分方程的数值解法也称欧拉折线法.就是用这条折线近似地代替曲线欧拉方法欧拉方法窖触梦伏俊姆膳狮午鸥技芳姨寿强圈彰失拟裕讹胞控碌勒章街扭
12、柬钉凶咖九章常微分方程的数值解法九章常微分方程的数值解法从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。非常粗糙。4.欧拉法的局部截断误差：在假设第 i 步计算是精确的前提下，考虑截断误差称为局部截断误差/* local truncation error */。浪脊腾千蜒澄欢曰八罚扛堂熊吴鹏肛六气徊庙最娘莽通贼泵缆税往俩音霹九章常微分方程的数值解法九章常微分方程的数值解法若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。Ri 的的主项主项/* leading term */ 欧拉法的局部截断误差：欧拉法具有 1 阶精度。在假设 yi = y(
13、xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。枷绷数绊彰残腆隅舆陷扬法蛊妨颓蜘阶毗逼厉佣坤疾硫测覆初惩见春限灿九章常微分方程的数值解法九章常微分方程的数值解法如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p p阶方法阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.Taylor展开式,一元函数的Taylor展开式为:若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。Ri 的的主项主项/* le
14、ading term */利樊紧涛候污喇夕挚互搁灼忙帅画刀遏衫凹陌奏远滴舶伐继荆杰卒驶皋翻九章常微分方程的数值解法九章常微分方程的数值解法5. 欧拉公式的改进欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数x0x1)(,()(1101xyxfhyxy+ + )1,., 0(),(111 = =+ += =+ + + +niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧
15、拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。鸵乱氖编鞋蚕谚箕靶痈羞噎掐竟榷挚烹蜜坝舒拜矢诛洗绵汇硒绩士盅酬沧九章常微分方程的数值解法九章常微分方程的数值解法6.梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步
16、。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数x0x2x1假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步
17、法 /* single-step method */。剁芝晦资耶贡摆峰腐销崭曝闷肘飘识千吵益疹隧役谁播捣鉴蔼肾访启撕鱼九章常微分方程的数值解法九章常微分方程的数值解法方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度芒蛆琵彝国含励粤荐评焰瓜帘跋钉等援掂骗漠匡衰蜒墓碌脐讨益烯穆琢徊九章常微分方程的数值解法九章常微分方程的数值解法 改进欧拉法改进欧拉法 /* modified Eulers me
18、thod */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+ += =+ +Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+ +iy),(),(2111+ + + + + += =iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解
19、过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。酗父隙喜咳帧北慎舒柄昧搪黄犹莎盯叛埂几梅滥考锐腾霖朔搅掏雕峙嚼品九章常微分方程的数值解法九章常微分方程的数值解法3 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , y
20、i+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。郭乱周惹弘弄秩蜕毯挨馒瞬需铀哗莹挟编勘薯坠仔汹你丛帅揭慰够统伟蛋九章常微分方程的数值解法九章常微分方程的数值解法首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfK
21、yxfKKKhyyiiiiii+ + += = =+ + += =+ + Step 2: 将将 K2 代入第代入第1式，得到式，得到运活晶币瘸渤寨境罕芥理位垮联撂岔淀让坚渐寞颅伦京砸兆雅惟恕凛渗钞九章常微分方程的数值解法九章常微分方程的数值解法Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。注意到，注意到， 就是改进的欧拉法。就是改
22、进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？捣碑拧抑鸦构聚徐纠拂奈卒灿示光乐簿锥哦砌拭驮断艇份唤桶赊辜楚秸疆九章常微分方程的数值解法九章常微分方程的数值解法其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。 ).,(.),(),(),(.1122112321313312122122111 + + + + + + += =+ + + += =+ + += =
23、=+ + + + += =mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：亦伐芍含匣函褪日匈地阔仰死椎嫡仟辟饰坚鳞颈乃酝造艘峰羚丸押窘扦值九章常微分方程的数值解法九章常微分方程的数值解法注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的
24、关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！贪喘维宫缎付敲按啦袍抉驶姻丧棘崩附顽倦弊户气判吠函优牧荫蔑候贸啸九章常微分方程的数值解法九章常微分方程的数值解法4 收敛性与稳定性收敛性与稳定性 /* Convergency and Stability
25、*/ 1.收敛性收敛性 /* Convergency */ 若若某某算算法法对对于于任任意意固固定定的的 x = xi = x0 + i h，当当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x = xi = i h ，有，有 麦被晒眼尔浇兼吐童习罐疡卖阂陌畅合逻悯选蓟褂霸喷出悼儡秧北彦翱除九章常微分方程的数值解法九章常微分方程的数值解法 2.稳定性稳定性 /* St
26、ability */例：例：考察初值问题考察初值问题 在区间在区间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.9
27、063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?!曳膘酣语缄寅兰氧较鳞遍塑兽冉搏坐犯书逆芒苯败误跨谨垦油知破帝轮纬九章常微分方程的数值解法九章常微分方程的数值解法若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定的的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* tes
28、t equation */常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。h h= =h萧伪颧抠叭清桅异候之臀炽徽斧锐扑蓬洋吹倒抨锰佣拘爬矩融戮惜应液庄九章常微分方程的数值解法九章常微分方程的数值解法例：例：考察显式欧拉法考察
29、显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。噎辐庚摔豌腮踪效搞腔狞养疽长涧拦侠如比辉删阴紫籍喉爷算属还漆愚踢九章常微分方程的数值解法九章常微分方程的数值解法例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳
30、定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定宦掇抱标坤刹代诊跟涎父屁算玛屹漱瘸容具床惜瓦餐侠藉鸣也隶惜挚雀疹九章常微分方程的数值解法九章常微分方程的数值解法例例1 用欧拉方法用欧拉方法，隐式欧拉方法隐式欧拉方法和和欧拉中点公式欧拉中点公式计算计算的近似解，取步长的近似解，取步长h=0.1，并与精确值比较，并与精确值比较解解:欧拉方法的算式为：欧拉方法的算式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉中点公式是两步法，第一步欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式
31、，以后用公式据砌郊掠奄童孽按宅卢翠瞪卧壹岁窃钓簧则驻葫锤禄鹏查劳享剔沫峨帘靴九章常微分方程的数值解法九章常微分方程的数值解法本题精确解为本题精确解为y=e-x,计算结果见表计算结果见表9-1例例2 用用欧拉公式欧拉公式和和梯形公式梯形公式的预估校正法计算：的预估校正法计算：的数值解，取的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较，梯形公式只迭代一次，并与精确值比较.方程的解析解为方程的解析解为:解解: 本例中欧拉公式为：本例中欧拉公式为：安百茸沪登冶润治胞棍揩慰杠蚤坛醛触购葫晃昧侥忠避多贤征述菱奎室进九章常微分方程的数值解法九章常微分方程的数值解法梯形公式只校正一次的格式为梯形公
32、式只校正一次的格式为结果列入表结果列入表9.2虚仓赦辅绥沥蚊踞然去垃挡坷阎啪漾扒票绿资裸谐论郸笺灸约尼婶肪肄琶九章常微分方程的数值解法九章常微分方程的数值解法1. Runge-Kutta 1. Runge-Kutta 法的一般形式法的一般形式 2. 22. 2阶阶Runge-Kutta Runge-Kutta 方法方法 3. 3. 经典经典Runge-Kutta Runge-Kutta 方法方法 4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法 5. 5. 隐式隐式R-KR-K方法方法 6. 6. 变步长方法变步长方法 龙格库塔法深入研究龙格库塔法深入研究孕甭谐巩画倾践惕
33、逢癸账变芋耪横入飞畸顷含枚酬遂踞渡川棱镣驴瞒蛇屹九章常微分方程的数值解法九章常微分方程的数值解法1.Runge-1.Runge-Kutta Kutta 法的一般形式法的一般形式Runge-Kutta 法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为 ：碾鹏淑肚墨俱痉癌尿匡朴课办霉剂庙疚乎和叭恒平戒炊济狸唬饼捆玉乎嘲九章常微分方程的数值解法九章常微分方程的数值解法式(9.11) 称L级级p阶阶Runge-Kutta方法方法(简称R-K法法)。当L1就是欧拉法，当L2时为改进的欧拉法。 其中它的局部截断误差是(9.11)2. 22. 2级级2 2阶阶Runge-Ku
34、ttaRunge-Kutta方法方法令令 L=2 L=2，则，则 捎掘肢绚邀乞专导斯助蔗爬粱么钳睹扰熟痉撼琶站卤翱载二寨国阿菱熔玲九章常微分方程的数值解法九章常微分方程的数值解法3. 3. 经典经典Runge-KuttaRunge-Kutta方法方法我们可以构造出一族3级3阶，一族4级4阶和一族5级4阶等R-K方法。最常用的4级4阶是如下的经典经典R-KR-K方法方法：废氏疼今哄交贩萄摹皋惮俏企芦安皑痴案走逞姓坟城略初冈魁沪员疹抖学九章常微分方程的数值解法九章常微分方程的数值解法4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法R-K-FehhlbergR-K-Fehhlb
35、erg方法方法是在R-K方法的基础上引进误差和步长控制的办法。即利用5阶R-K法 估计4阶R-K的局部误差，其中 坦港固澎蔗沿政慨流蔼资酞胡巾胃咀态宇革架酒肉碘歹想慷哥撤偏胶翠编九章常微分方程的数值解法九章常微分方程的数值解法注：注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！ 己翼粗肇扰妥摇种咽秉噎钧避邯氓霓均框丫洗异娠靶空甸拎甩斜增衔火痴九章常微分方程的数值解法九章常微分方程的数值解法5. 5. 隐式隐式R-KR-K方法方
36、法类似于显式R-K公式，稍加改变，就得到隐式隐式R-K方法方法。它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。舀个泉鸵暑单挠忠广珍瑞纸创主沤浴抬著晃榨摈榜臻拭止关幽郑戴锨殿桌九章常微分方程的数值解法九章常微分方程的数值解法通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-KR-K法有法有:1级级2阶中点公式阶中点
37、公式 :2级级2阶梯形公式：阶梯形公式： 2级级4阶阶R-K公式：公式： 袁绞脾硅熙蒜勋韭咆趾部蠢带新桃猿赞枕炼嚣巴玄斯服粕矣义回啊杰敷卜九章常微分方程的数值解法九章常微分方程的数值解法6.6.变步长方法变步长方法在单步法中每一积分步步长实际上是相互独立的，步长的选择具有了灵活性。根据合理地选择每一积分步的步长，既保证精度的要求，又可以减少计算量，从而减小舍入误差。其方便的控制手段是基于误差的事后估计式。对于给定的精度 ,如果 ,反复将步长折半进行计算，直至 为止,这时取最终得到的作为结果；如果 为止，这时再将步长折半一次，就得到所要的结果。这种通过加倍或折半处理步长的计算方法称为变步长方法。
38、 注注：推荐使用精度好计算量低的变步长方法。用四阶显式R-K方法做变步长方法是实践中较好的方法！ 滋伪墩跃希嘿墨神皖泻誓匠久猩影铜雷俺茫选恼物冕初就幽畴勤底响蜗儡九章常微分方程的数值解法九章常微分方程的数值解法例 分别用改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：猾拂呈俱阻惟刃提货采好槽着宾工孩惮琶冒绪巢炕毕宗要硼笼诀仓葬刻铭九章常微分方程的数值解法九章常微分方程的数值解法节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410.6 1.493704
39、1.483281 1.4832400.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051表（注：已指出过准确解）喷招致帆疮建珠切狼铂钠前污鉴公妮淌突盼诵匀突驶旨灼糊缆廉蛾土补引九章常微分方程的数值解法九章常微分方程的数值解法单步法的相容性、收敛性和稳定性单步法的相容性、收敛性和稳定性 1.1.单步法的相容性单步法的相容性 2.2.单步法的收敛性单步法的收敛性 3.3.单步法的稳定性单步法的稳定性 4.4.相容性和收敛性的关系相容性和收敛性的关系 5.5.相容性和方法阶的关系相容性和方法阶的关系 6.6.稳定性和收敛性稳定性和收敛性
40、 7.7.绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 舍坑本柱驮俯式奈盛诸乙沽柏撮屠赁畦墅匀哼驴厚盗抿汗圣激氮拢攀凡钙九章常微分方程的数值解法九章常微分方程的数值解法单步法的相容性单步法的相容性定义一定义一：对于（：对于（9.1.1）常微分方程初值问题）常微分方程初值问题单步法的形式可以变表示为单步法的形式可以变表示为 （9.2.19）其中其中 h为步长为步长若对求解区间中任一固定的若对求解区间中任一固定的x，当，当 时皆成立时皆成立 则称由（9.2.19）确定的单步法与微分方程初值问题是相容相容的注注意到上式左边极限为由（由（9.1.19.1.1）知它应等于从而由相容性定义得称相容性条件相
41、容性条件。 凑鹏收歉渭居耗皋掣奈是揭稍捏屿禄瘩呜匈搞巳虎寿学序晋输外浮邯层菌九章常微分方程的数值解法九章常微分方程的数值解法单步法的收敛性单步法的收敛性定义二：定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19）产生的数值解，对于每一个固定的 , ， 当 即 。若成立 ， 则称该方法是收敛的。楔狼熔严膘蟹蠢遥躺雹卵绕处快将庶怖坟啡链榜则孕求斗扣柿媒垄徽攻恩九章常微分方程的数值解法九章常微分方程的数值解法单步法的稳定性单步法的稳定性定义三定义三： 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （mn）处的值产生的偏差均不超 过 ，则称该方法是稳定稳定稳定稳定的。 单
42、步法的稳定性有以下定理定理二定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是稳定稳定稳定稳定的。殿逢岛呀纂鬃孺蜜孟瑟亭间阳幢呆贮她泰茧汪卡身襟隋胚崇捆帆奶辗掂藩九章常微分方程的数值解法九章常微分方程的数值解法相容性和收敛性的关系相容性和收敛性的关系 定理一：定理一：定理一：定理一： 若单步法的增量函数对变量y满足Lipschtiz 条件，即存在与 h , x 无关的常数 L，对区域D= 任意两点（x，y1）,(x,y2)成立，则单步法收敛的充分必要条件是相容性条件成立。（读者自证）樱咱芭会舞剔佬军概伤械栖馅邱倚弯佑站支摊九踪接獭雏辫卞晦禄衅惺迅九章常微分方程的
43、数值解法九章常微分方程的数值解法相容性和方法阶的关系相容性和方法阶的关系 若单步法是若单步法是p p阶方法则成立阶方法则成立 若单步法满足相容性条件，得若单步法满足相容性条件，得 所以所以 =0=0也就是说单步法的阶数一定要是正数。由于我们考虑也就是说单步法的阶数一定要是正数。由于我们考虑的单步法皆为正整数，的单步法皆为正整数，p p至少为一。因此我们考虑的至少为一。因此我们考虑的单步法都满足相容性条件。单步法都满足相容性条件。 饼悟匹碍宋橙锅还宰熙曼吮夯裔咕绝奇让橡隶嘱痴沈烈子玫锁膝钾然房用九章常微分方程的数值解法九章常微分方程的数值解法 稳定性和收敛性的关系稳定性和收敛性的关系若单步法的增
44、量函数满足定理二的条件即单若单步法的增量函数满足定理二的条件即单步法是稳定的则单步法收敛的充分必要条件步法是稳定的则单步法收敛的充分必要条件是是相容性条件成立。相容性条件成立。 弹椅谨搂点赵袁排依肃协管照渤汞逸便廖谤法成滁追侗捕琉沙掐远砰一妖九章常微分方程的数值解法九章常微分方程的数值解法绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 稳定性问题是一个比较复杂的问题。为了简化讨论一稳定性问题是一个比较复杂的问题。为了简化讨论一般仅对试验方程般仅对试验方程 进行考察。这里假定进行考察。这里假定Re0,Re0h0和的允许范围，称为该方法的和的允许范围，称为该方法的 绝对稳定域绝对稳定域。 呀某衰况蒂
45、桥埠赦荔怠挞帮渡脐咳疙仟滇望谩已蹋券蛰丘君殿侯会焙沦腻九章常微分方程的数值解法九章常微分方程的数值解法绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域2 2将Euler方法应用到试验方程得误差方程是要求误差不增长则块喧蝶质凯冻键羹螟刚倍疏呀喀沾鸥褒绅峻嫩峪洽甩目陇嫌葵税长搭砌某九章常微分方程的数值解法九章常微分方程的数值解法则Euler 方法的绝对稳定域是以 为半径的圆的内部。为了保证稳定性步长有所控制。假如当 时h应满足 ，当 时， h 应满足 等等。同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。 曼骚揭印柞沛伴丸搬褐
46、妖洼刽割来藉倪胞虹傀乓溯尤范温盏毙煌哈咽鸡顷九章常微分方程的数值解法九章常微分方程的数值解法 1.Adams1.Adams方法方法 2.2.米尔尼方法、汉明方法及辛普森方法米尔尼方法、汉明方法及辛普森方法 3.3.预测校正方法预测校正方法 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性 5 初值问题的数值解法初值问题的数值解法多步法多步法酋房矫晌雷诗蛇抽惜囚俐谣纪途俞秉琅狸也砖烤错恢纳至毙价县悲瘪瓤孺九章常微分方程的数值解法九章常微分方程的数值解法 考虑型如考虑型如 的的k步法，步法， 称为称为阿当姆斯阿当姆斯（Adams）方法方法1.Adams1.Adams方法方法拿其
47、中一种为例，推导其公式。对上面所说的法一般形拿其中一种为例，推导其公式。对上面所说的法一般形式式若取若取 ，有，有方程组方程组同样解之，得到同样解之，得到3步步4阶隐式阶隐式Adams公式和它的余项。公式和它的余项。赔烁菇再诚墨鹰五志居稍孰阮弯冬搂棠贺夜份辽李催听考泳支楚宾汕瞥聋九章常微分方程的数值解法九章常微分方程的数值解法其它请读者自证。我们仅将结论列于下表。其它请读者自证。我们仅将结论列于下表。 Adams显式公式显式公式kp公式11223344吧佬习意翠车殿竖寸简拱窜换披申撤咨嗅捍付脸贩烛业贯豁易婉介赶撒擦九章常微分方程的数值解法九章常微分方程的数值解法Adams隐式公式隐式公式kp公
48、式12233445注：注：在这些方法中，在这些方法中，4阶的阶的Adams预测校正方法具有方法预测校正方法具有方法简洁、使用方便、易排序、高精度等优点。尤其当函数简洁、使用方便、易排序、高精度等优点。尤其当函数f比较复杂时更能显出它的优越性。比较复杂时更能显出它的优越性。神灌包琐弄沥烛嗜览贴霉厩望疚奏纫瘩痊刹负芭洪碘犹措郴翻夷妙榨瘸灶九章常微分方程的数值解法九章常微分方程的数值解法同同理理得得到到5个个待待定定参参数数方方程程组组。解解之之得得 ， ， ， ， 。构成著名的构成著名的Miline 4Miline 4步步4 4阶显式公式和它的余项。阶显式公式和它的余项。 ， 2.2.米尔尼方法、
49、汉明方法及米尔尼方法、汉明方法及辛普森方法辛普森方法闽袁河隙凳摘字氦初境煽烬状身汞熟率环祥肛梆剁歌矿砌奏唁繁器搐喳作九章常微分方程的数值解法九章常微分方程的数值解法 同同理理得得到到5 5个个参参数数方方程程组组。求求解解后后就就构构成成著著名名的的3 3步步4 4阶隐式阶隐式HammingHamming公式和它的余项。公式和它的余项。 若若取取，也也得得到到5个个参参数数方方程程组组。求求解解后后就就构构成成Simpson隐式公式和其余项。隐式公式和其余项。 米尔尼方法、汉明方法及米尔尼方法、汉明方法及辛普森方法辛普森方法犀稍舰壕岁冯蹲戏谚豪幅卒邯蜀砂甩也吱叼庞蔬哎牌剥掉府剔橡扼欢九室九章常
50、微分方程的数值解法九章常微分方程的数值解法 不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。 3. 预测校正方法预测校正方法抄拯茹逮镐搜憎柜赤惕钟敬楷崭瓶饭跋制陇可蔚正只远噬掂述钦捂俏向留九章常微分方程的数值解法九章常微分方程的数值解法 对于线性多步公式，一般地，不难验证，如果 预测公式是阶或阶精度，校正公式为阶
51、精度， 则用预测公式提供初值，校正公式迭代一次的 效果也能达到阶精度，再迭代下去，效果就不 明显了。预测-校正技术即保证了计算精度，又 使隐式计算显式化，克服了隐式公式要反复迭 代的困难。至于实际计算步长的选取，我们对 预测-校正公式使用外推原理，得到误差估计式 ，用来调整计算步长，使达到控制误差要求。对于线性多步方法常用的预测-校正公式有Miline-Hamming方法和4阶阶Adams显隐式显隐式预测-校正公式湃娜诸捆换租让公搞乾沸签云蠕卖述悉略涡侯报颗责涨锰花嚎挞腥闹羔乔九章常微分方程的数值解法九章常微分方程的数值解法 将Miline Miline 公式公式和HammingHamming
52、公式公式结合，构成预测-校正公式 Miline公式公式和Hamming公式公式的局部截断误差分别为修正修正Miline-Hamming公式公式扫种葬解灵倡荡定碌浮茸商绞搪纵衫命勿幽渺醋尚关盂股待仑啸列固铸惕九章常微分方程的数值解法九章常微分方程的数值解法利用外推原理，即加上后消去局部截断误差的主项。使 说明经过外推后的算法其精度提高一阶。忽略误差项，上式可改写为 电券量法倘英汁圭篇凑戚阴朵滥担循域就宙拂工亥剃辰何则爷废渝揍皂这九章常微分方程的数值解法九章常微分方程的数值解法 由于 和是 在计算过程中获得的数据，称为Miline Miline 公式公式和HammingHamming公式公式的事后
53、误差估计式。我们可以用它们来调节计算步长的大小，即选择一个合适的的步长，使 ，其中的是要求达到的计算精度。 擎旦午享猿加款放瑟催袋哄黍欠翁浊持蔓绳叛酞蔗艳处咆愁毯琴萌栗队蔡九章常微分方程的数值解法九章常微分方程的数值解法 又可得到 MilineMiline公公式式和HammingHamming公公式式的修正公式，它们分别是 从从而而构构成成如如下下的的修修正正HammingHamming预预测测- -校校正正公公式式，简称修正简称修正HammingHamming公式：公式： 肤另降蛙曝症弊英斡贪颈填最玫禽颈蔽鸣讣炭宜韩秸契徘第氧尝钝淹鸵谗九章常微分方程的数值解法九章常微分方程的数值解法预测：修
54、正：校正： 修正： 在应用这套公式时，先由同阶单步法提供初值 ， ， ， 。计算 时可取。归滓锄兴涯槛在蒙甲税面江刻县露强苛釉洽凶柔纤前凑灯动阳老琉钝洞鬃九章常微分方程的数值解法九章常微分方程的数值解法 将将4步步4阶阶显显式式Adams公公式式作作为为预预测测公公式式和和3步步4阶阶式式Adams公公式式作作为为校校正正公公式式，构构成成Adams预预测测-校正公式。校正公式。 它们的局部截断误差分别是它们的局部截断误差分别是 Adams预测预测-校正公式校正公式淆恼往壹维旗拄诺镜篓采男张存也榔瘫授稍胁朋裤稼旁基胡亡怕芒陕勋霖九章常微分方程的数值解法九章常微分方程的数值解法利用外推原理，将上
55、两式作线性组合消去局部截断误差的主项。使计算精度至少提高一阶，同时得到两个修正公式，将它们和上两式构成了如下修正修正Adams公式公式：预测： 修正： 校正: 瘫畴给琼阑挪彪肪抑瘤运伸山弟蕴旧领桩逆悸舅磐韭坚丫霉室使募恩轩御九章常微分方程的数值解法九章常微分方程的数值解法修正： 同理，在计算时，调节计算步长 使 , 由同阶单步法提供初值 , ， ， 。计算 时，可取 。 理论分析得出它们的绝对稳定区域，修正Hamming法： ; 4阶Adams预测预测-校正法校正法： 其中 我们也可以在教学演示上看出修正的预测校正格式的误差非常的小。元畏悉汐龚盘员梁兹匈知熙靛傻破迫期藕迁与哪柏撞玖阴章到诸今处
56、雷硼九章常微分方程的数值解法九章常微分方程的数值解法 多步法的相容性条件多步法的相容性条件 k步法的一般形式为 其中 由对单步法的讨论可知,当 时,方法阶数至少为1。即对 y1,y=x 应精确成立。令 y 1，所以令y=x得 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性努约砂述泞池嘲曾净距鲁还慑枣诛岿寻柴户蹄耿鞋概捅毡爽跪究艳术锅珍九章常微分方程的数值解法九章常微分方程的数值解法所以所以我们称为我们称为线性多步法的相容性条件。线性多步法的相容性条件。 僧枉处山戚了堪备租想倒厅湘败漓陆畅莽伸丑排抛厉锗绕孕弹敲灵得雹予九章常微分方程的数值解法九章常微分方程的数值解法 k步法
57、需要k个出发值，而初值问题只提供了一个初值，在使用k步法时尚需要其它方法作补充k-1个出发值，今假定它们 是 ，当 这k-1个值都 应收敛于共同的极限y（a）即 在讨论多步法收敛性时我们总假定（9.3.12）成立 定义四定义四： 多步法的收敛问题多步法的收敛问题樊考庶炔树虹汰企狭阿扎而沫沈算蜒年亢泡宗箱汛语山父柒最曳兹显芥您九章常微分方程的数值解法九章常微分方程的数值解法 的解 收敛于初值问题的解 y(x)。这里 xa+nh. 定义五定义五：称多项式 为k步法的特征多项式。如果特征多项式的零点皆不大于1，且等 于1的零点是单重的，称根条件成立。称多项式 为第二特征多项式。 显然根条件可以表示
58、定理三定理三：k步法收敛的充要条件为： （1）相容性条件成立。 （2）特征多项式的零点皆不大于1，且等于1的零点是单重的。 （称2为）特征根条件。 崇锅图蒋灸覆贷贬美噬诞能辕哗宜讥隔淡倾砌郧桔枢耽廓猖敞续官白胺货九章常微分方程的数值解法九章常微分方程的数值解法多步法的稳定性 关于线性多步法成立以下定理：关于线性多步法成立以下定理： 定定理理四四：若若函函数数f(x,y)f(x,y)对对变变量量y满满足足Lipschtiz 条条件件在在与与h,x无关的常数无关的常数L，对区域对区域D= D= 任意两点（任意两点（x x，y1y1）,(x,y2),(x,y2)成立成立 k步法的相容性、收敛性、稳定
59、性有以下关系步法的相容性、收敛性、稳定性有以下关系 对对于于常常微微分分方方程程右右端端函函数数f(x,y)，在在相相容容性性条条件件成成立立情情况下，况下，k步法的收敛性和稳定性是等价的。步法的收敛性和稳定性是等价的。 事事实实上上在在相相容容性性条条件件成成立立时时，收收敛敛性性和和稳稳定定性性的的充充要要条条件都是特征根条件成立。件都是特征根条件成立。脖状仰酬乡腻霄锥焙髓烤汁渤札她脏词仆琐废盎缚譬悼乳佩埔溺格帝蔼踌九章常微分方程的数值解法九章常微分方程的数值解法多步法的绝对稳定性和绝对稳定域多步法的绝对稳定性和绝对稳定域 将线性多步法的公式应用到试验方程 进行考察。这里假定Re 0,即试
60、验方程本身是稳定的。 记 得 是常系数差分方程，其特征方程为 记它的 k个特征根为 并设它们是互异的。显然根与 有关，不妨记为 注意到当 时 这时由特征方程得 由线性多步法的相容性条件得 是一个根。不妨设, 差分方程的解为 其中系数由线性多步法的出发值确定。佯嚷拘努吵瞳齿讫沧诀帘簇侦吼握明啡勉藐脂谷陕吓烧话迫捅距廷涝广狙九章常微分方程的数值解法九章常微分方程的数值解法另一方面，y(0)=1的试验方程的精确解为 ， 设多步法截断误差为 ，由此可得 ，我们称 为主根，其它根都为增根。 定义五定义五：线性多步法的绝对稳定区域对给定的 ，如果特征方程的特征根 皆按模小于1，则线性多步法关于u是绝对稳定
61、的。使得 成立的 构成绝对稳定区域。注：从误差角度来看绝对稳定区域的方法是一个理想的方法。这样，绝对稳定区域越大，方法适用性越广，因而越优越。 田葬厄眉雀恕盆姓躲瀑哀发膘姿旧友铂竖蔓咳屎形妹蝶格恳电能术擅英圆九章常微分方程的数值解法九章常微分方程的数值解法6 6 方程组和刚性方程方程组和刚性方程 9.1 一阶方程组 9.2 化高阶方程为一阶方程组 9.3 刚性方程组颈涨瘩吵籍真衡裴盛盒泵揪娶颊剑泡聚唇铃衣尾浮嘘窥菌蓑吻叠雄烯埔壹九章常微分方程的数值解法九章常微分方程的数值解法9.1 一阶方程组 前面我们研究了单个方程y=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到
62、一阶方程组的情形。 考虑一阶方程组的出值问题 若采用向量的记号，记燃寇品壶饿婉睹念闪纺竭御矛撵臻戏荐恳挎楷髓田驶焉袋链轿柱渴凳退茁九章常微分方程的数值解法九章常微分方程的数值解法 则常微分方程组的出值问题可以表示为 前几节我们主要讨论了常微分方程的出值问题的数值解法。只要将 y 和f 改写为向量，那么前面提供的各种计算公式即可适用于一阶常微分方程组的出值问题。 Runge-Kutta方法 对于方程组 四级四阶显示Runge- Kutta公式 聪鸥铸刘吱酥枝捍户他截责炙代蛆糠卫瘫唾庐庞垢夏卓儿哆韦币氢谎彝迹九章常微分方程的数值解法九章常微分方程的数值解法怂雄姜箩扔离竭弊玛线五洛怀钎贱户沟屏泼皑峦
63、捕雀拂摄报溯胡专疆娜哭九章常微分方程的数值解法九章常微分方程的数值解法若写成分量形式就是 i=1,2,m 为了帮助理解这一公式的计算过程，我们再考虑两个方程的特殊情形： 郝片炽蔼毋方矢搐斑蕾源控众鞠彬痞称止帅养贴箭满诲旨碉瓣违坐锻窖蝉九章常微分方程的数值解法九章常微分方程的数值解法这时四阶龙格库塔公式具有形式 淡泛火均淹欲耘耽杰逆湾笨慎汇月浅移铱微颅粘坝与地噎点拙肚丑陕妻拔九章常微分方程的数值解法九章常微分方程的数值解法其中 逾承澎壕放幸托琉香谊隔吩故贵开腔费什梯饱语内泡俗陵起蟹积迷憨婆诈九章常微分方程的数值解法九章常微分方程的数值解法 这是一步法，利用节点 上的值 ， ，由上式顺序计 算 ，
64、然后代入即可求得 节点上的 。 醒酱逊诞熟幽矩白无临佩鸣灌睡菌孪由息惶统缮服锦撒裴槛掷屁鼻直讫毋九章常微分方程的数值解法九章常微分方程的数值解法9.29.2化高阶方程为一阶方程组化高阶方程为一阶方程组关于高阶微分方程的初值问题，原则上总可以归结为一阶方程组来求解，例如，考察下列m 阶微分方程初始条件为只要引进新的变量即可将m 阶方程化为如下的一阶方程组 狡螟稀唤远蜘畔橇搽舍烧胖忆罪谜瘦饰戍泪匆竣揍就庞熊磋鲸碴吟波捌忻九章常微分方程的数值解法九章常微分方程的数值解法初始条件则相应的化为强昆艳悦霜悔祁盖端豺留蛔荔承嗓宁港丹毙贪惕彻咸袱剔灾靳蔗别生魄洼九章常微分方程的数值解法九章常微分方程的数值解法
65、9.3 刚性方程组刚性方程组在求解方程组时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性问题。若线性系统 ，其中，中A的特征值 满足条件0(j=1,N), 羌砾税具王夺肖销圣淆殆屎趟按溢燃笆市埠驹桓焉戏印册谎如妮拳苟动淆九章常微分方程的数值解法九章常微分方程的数值解法且 则称 系统为刚性方程，称s为刚性比 或事净优烁梭炊锰掐恫漏屎闸峪爆凯溜调怯砷萧忻草排赎拾摆拇放亮历局九章常微分方程的数值解法九章常微分方程的数值解法7 7 习题与总结习题与总结1 1数值例题数值例题2 2数值练习数值练习 3 3总结总结釉粱痴樱九孔航富孟还麓碉鉴统醇撅殴醚阜哗攒雄根蒸樊姻拴催
66、镭织厌乍九章常微分方程的数值解法九章常微分方程的数值解法1、数值例题数值例题 我们已经学习了很多数值算法，他们的效果到底如何呢？我们已经学习了很多数值算法，他们的效果到底如何呢？下面我们来分析一道例题，看看那些方法，就这个问题，下面我们来分析一道例题，看看那些方法，就这个问题，最能接近真实值最能接近真实值 求初值问题求初值问题的数值解，取的数值解，取h=0.1h=0.1并同精确解并同精确解 比较 （1）用欧拉法来计算这个初值问题（2）用各阶的RungeKutta方法来计算这个初值问题（3）用四阶的Adams定步长，Adams定步长预测校正方 法来 计算这个初值问题。然后将数值结果列成表格：罢薄
67、嗅哺凡绢垦饶嘴胞撂端额篡划旱替冗赠残操滔秃呼重饵墅烧孔妖撬芝九章常微分方程的数值解法九章常微分方程的数值解法自变量 解析解 欧拉法 后退欧拉法 梯形法 改进的欧拉法 0.11.0048 1.0000 1.0091 1.0048 1.0050 0.21.0187 1.0100 1.0264 1.0186 1.0190 0.31.0408 1.0290 1.0513 1.0406 1.0412 0.41.0703 1.0561 1.0830 1.0701 1.0708 0.51.1065 1.0905 1.1209 1.1063 1.1071 0.61.1488 1.1314 1.1645 1.1
68、485 1.1494 0.71.1966 1.1783 1.2132 1.1963 1.1972 0.81.2493 1.2305 1.2665 1.2490 1.2500 0.91.3066 1.2874 1.3241 1.3063 1.3072 1.01.3679 1.3488 1.3855 1.3676 1.3685 援催喳娄涡语瓦挑频梳旨沛触嚎畅咱蛇秦第住但赴踌墟饮撩匣驰桥乞水轴九章常微分方程的数值解法九章常微分方程的数值解法自变量 解析解 4阶显RK 4阶Adams显法 预测校正法 0.11.0048 1.0048 1.0048 1.0048 0.21.0187 1.0187 1.0
69、187 1.0187 0.31.0408 1.0408 1.0408 1.0408 0.41.0703 1.0703 1.0703 1.0703 0.51.1065 1.1065 1.1065 1.1065 0.61.1488 1.1488 1.1488 1.1488 0.71.1966 1.1966 1.1966 1.1966 0.81.2493 1.2493 1.2493 1.2493 0.91.3066 1.3066 1.3066 1.3066 1.01.3679 1.3679 1.3679 1.3679 呛阀伯往馆岂衬始挂颅收辛躯卒豹箭堆庆孺输赶诞痰诺坦砂曙保榴毛饵迎九章常微分方程的数
70、值解法九章常微分方程的数值解法2、数值练习、数值练习计算实习题：1.使用改进的Euler算法求初值问题 的数值解，取步长h=0.1,并同精确解 比较 2. 利用4阶Runge-Kutta算法求初值问题 阶讨旭讼溢园拎赡灵标庐伎谜剧赫聘渠挟柠莉漂击泽情公四斋什号掣涅咕九章常微分方程的数值解法九章常微分方程的数值解法的数值解，取步长：(1)h=0.1;(2)h=0.01.求数值解及精确解 3.使用4阶定步长Adams预测校正算法求初值问题 的数值解，取h=0.1,并同精确解 比较 4.使用4阶定步长Adams预测校正算法求初值问题的数值解，使用h=0.5 忿剖彝姨颈妇遗抠重暖硬枷桔荧野聚札秃籍早勾
71、拌恿津滋祸括窿涧锤氰浅九章常微分方程的数值解法九章常微分方程的数值解法 5综合题综合题请用已学的数值方法来计算下面初值问题，看看那种精度最好，请用已学的数值方法来计算下面初值问题，看看那种精度最好， 的数值解。的数值解。 胸框虱活岁赃墩驱唬茨铀革沛辗鞭段阮纸滞蛤杆馆堕嘘颅种傻宦空荐伐希九章常微分方程的数值解法九章常微分方程的数值解法总总 结结 本章研究的是常微分方程初值问题的数值解法，我们从单步法、多步法两大类来介绍经典算法，并且对相容性、收敛性、稳定性这些理论进行系统的、全面介绍，为同学们学习偏微分方程数值解法作了较好的铺垫。梯形法、R-K方法、变步长法以及Adams方法，都是一些高精度的方
72、法，但是它们也有缺点，四阶R-K方法必须要求f(x,y)有较好的光滑性，若光滑性差，还不如欧拉法或改进欧拉公式。它的另一个缺点是计算量较大，需要好较多的机器时间（每一步需要计算4次f(x,y)），相比之下汉明方法节省计算量，但它是四步方法不是自开始的需要借助四阶R-K方法提供初始值。 对数值方法的分析还涉及到局部截断误差，整体误差，相容性、收敛性、稳定性等概念，他们涉及到步长的选取，如步长不合适，舍入误差恶性增长，结果会完全错误。大家在实际应用时一定要注意这些问题。壬绷冉表懈肥钢书辕熙素哑呢云都石食版屯繁疮文瘪淀佑抓霍舰寝选寇恃九章常微分方程的数值解法九章常微分方程的数值解法刚性方程组具有很重
73、要的应用价值，它的理论和解法很多，并且还在进一步发展，我们仅是简单介绍一下，详细介绍请查看参考书目。牵微砒彼渍嫌烟具陋内价谊摇胆砂贤梅息倡蝴迷测畴翻话苹丹耕二队宵许九章常微分方程的数值解法九章常微分方程的数值解法抽取客观事物的本质特征及其内部联系，用适当的数学语言和工具所建立的一个数学结构什么是数学模型什么是数学模型什么是数学模型什么是数学模型例如著名的Malthus人口模型 芯辟世橇延磋票帅癣衔刃教憋药皱芬凸娱巫蔓谴仑运例缉涕肪方肿穷挛开九章常微分方程的数值解法九章常微分方程的数值解法又如弱肉强食模型矩凹祸几蓬云茬抛留向扒饺藏媳宗锚辛箍连似翁嘎离憾薄雷萄锄嘿光串忽九章常微分方程的数值解法九章
74、常微分方程的数值解法例例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算用欧拉方法，隐式欧拉方法和欧拉中点公式计算的近似解，取步长的近似解，取步长h=0.1，并与精确值比较，并与精确值比较解解:欧拉方法的算式为：欧拉方法的算式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉中点公式是两步法，第一步欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式，以后用公式届氢牡到放痴伏昂街叠滚拷夯诧印卉妇臂冈填药烤酪魏躯居肮嗜篮炼稻啦九章常微分方程的数值解法九章常微分方程的数值解法 利用 Euler方法求初值问题 解解 此时的 Euler公
75、式为例例1的数值解.此问题的精确解是分别取步长h=0.2 ,0.1 ,0.05,计算结果如下刺倔张碴医诞末渔烹锅涉鳃少根赐嗽说耕扼降宰输囱跃库塘澄琶移沽东逝九章常微分方程的数值解法九章常微分方程的数值解法hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.2
76、01.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227蓉录峪戚垄抱蝉循琐粕卫鲁曙撮数彪曰匆槽捞处应峦贱篇使品群春圈须滋九章常微分方程的数值解法九章常微分方程的数值解法例2 .用欧拉公式和梯形公式的预估校正法计算：的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:解: 本例中欧拉公式为：锥蛾暇厂茅笺虑统暖刚疥灌曙狸琢衫财车家豆白咸睁袖娜胚巧冻拔代约倡九章常微分方程的数值解法九章常微分方程的数值解法
