数值分析常微分方程数值解法省公共课一等奖全国赛课获奖课件.pptx

（Numerical Methods for Ordinary Differential Equations）第1页问题驱动：蝴蝶效应问题驱动：蝴蝶效应 洛伦兹吸引子洛伦兹吸引子(Lorenz attractor)是由是由MIT大学气象学家大学气象学家Edward Lorenz在在1963年给出，他给出第一个混沌现象年给出，他给出第一个混沌现象蝴蝴蝶效应。蝶效应。图图10.1.1蝴蝶效应示意图蝴蝶效应示意图 洛伦兹方程是大气流体动力学模型一个简化常微分方程组：洛伦兹方程是大气流体动力学模型一个简化常微分方程组：第2页 该方程组起源于模拟大气对流，该模型除了在天气预报中有显该方程组起源于模拟大气对流，该模型除了在天气预报中有显著应用之外，还能够用于研究空气污染和全球侯改变。洛伦著应用之外，还能够用于研究空气污染和全球侯改变。洛伦兹借助于这个模型，将大气流体运动强度兹借助于这个模型，将大气流体运动强度x与水平和垂直方与水平和垂直方向温度改变向温度改变y和和z联络了起来。参数联络了起来。参数 称为普兰特数，称为普兰特数，是规范是规范 化瑞利数，化瑞利数，和几何形状相关。洛伦兹方程是非线性方程组，和几何形状相关。洛伦兹方程是非线性方程组，无法求出解析解，必须使用数值方法求解上述微分方程组。洛无法求出解析解，必须使用数值方法求解上述微分方程组。洛伦兹用数值解绘制结果图伦兹用数值解绘制结果图10.1.1，并发觉了混沌现象。，并发觉了混沌现象。第3页1 1 引引 言言 微分方程数值解普通可分为：常微分方程数值解和偏微分微分方程数值解普通可分为：常微分方程数值解和偏微分方程数值解。自然界与工程技术中许多现象，其数学表示式方程数值解。自然界与工程技术中许多现象，其数学表示式可归结为常微分方程（组）定解问题。一些偏微分方程问题可归结为常微分方程（组）定解问题。一些偏微分方程问题也能够转化为常微分方程问题来（近似）求解。也能够转化为常微分方程问题来（近似）求解。Newton最早采最早采用数学方法研究二体问题，其中需要求解运动方程就是常微用数学方法研究二体问题，其中需要求解运动方程就是常微分方程。许多著名数学家，如分方程。许多著名数学家，如 Bernoulli（家族），（家族），Euler、Gauss、Lagrange和和Laplace等，都遵照历史传统，研究主要等，都遵照历史传统，研究主要力学问题数学模型，在这些问题中，许多是常微分方程力学问题数学模型，在这些问题中，许多是常微分方程求解。作为科学史上一段佳话，海王星发觉就是经过对常求解。作为科学史上一段佳话，海王星发觉就是经过对常微分方程近似计算得到。本章主要介绍常微分方程数值解微分方程近似计算得到。本章主要介绍常微分方程数值解若干方法。若干方法。第4页一一、初值问题数值解法初值问题数值解法1、一阶常微分方程初值问题普通形式、一阶常微分方程初值问题普通形式常微分方程数值解法分为常微分方程数值解法分为（1 1）初值问题数值解法）初值问题数值解法 （2 2）边值问题数值解法）边值问题数值解法第5页(2)普通结构方法：普通结构方法：离散点函数值集合离散点函数值集合离散点函数值集合离散点函数值集合+线性组合结构线性组合结构线性组合结构线性组合结构 近似公式近似公式近似公式近似公式2.迭代格式结构迭代格式结构(1)结结构构思思想想：将将连连续续微微分分方方程程及及初初值值条条件件离离散散为为线线性性方方程程组组加加以以求求解解。因因为为离离散散化化出出发发点点不不一一样样，产产生生出出各各种种不不一一样样数数值值方方法法。基基本本方方法法有有：有有限限差差分分法法（数数值值微微分分）、有有限限体体积积法法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。第6页(3)怎样确保迭代公式稳定性与收敛性怎样确保迭代公式稳定性与收敛性?3.微分方程数值解法需要处理主要问题微分方程数值解法需要处理主要问题(1)怎样将微分方程离散化，并建立求其怎样将微分方程离散化，并建立求其数值解迭代公式？数值解迭代公式？(2)怎样预计迭代公式局部截断误差与整体误差？怎样预计迭代公式局部截断误差与整体误差？第7页称称 在区域在区域D上对上对 满足满足Lipschitz条件条件是指是指:记记4、相关定义、相关定义第8页二、初值问题解存在唯一性二、初值问题解存在唯一性 考虑一阶常微分方程初值问题考虑一阶常微分方程初值问题/*Initial-Value Problem*/:则上述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上连续,且关于且关于 y 满足满足 Lipschitz 条件，条件，即存在与即存在与 无关常数无关常数 L 使使对任意定义在对任意定义在 上上 都成立，都成立，第9页 求函数求函数 y(x)在一系列节点在一系列节点 a=x0 x1 xn=b 处近似值处近似值 方法称为微分方程数值解法。方法称为微分方程数值解法。称节点间距称节点间距 为步长，为步长，通常采取通常采取等距节点等距节点，即取，即取 hi=h(常数常数)。称为微分方程称为微分方程数值解数值解。第10页三三、初值问题离散化方法初值问题离散化方法 离散化方法基本特点是依照某一递推公式，离散化方法基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右次序依次求出按节点从左至右次序依次求出 近似近似 假如计算假如计算 ，只用到前一步值，只用到前一步值 ，则称这类则称这类方法为方法为单步方法单步方法。假如计算假如计算 需用到前需用到前r步值步值 ,则称这类方法为则称这类方法为r步方法步方法。第11页2 2 欧拉方法欧拉方法 /*Eulers Method*/*Eulers Method*/欧拉公式欧拉公式(单步显示公式）：单步显示公式）：向前差商近似导数向前差商近似导数记为记为x0 x1第12页亦称为亦称为欧拉折线法欧拉折线法/*Eulers polygonal arc method*/在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是准准确确前前提提下下，考考虑虑截截断断误误差差 Ri=y(xi+1)yi+1 称称为为局局部部截截断断误误差差/*local truncation error*/。若若某某算算法法局局部部截截断断误误差差为为O(hp+1)，则则称称该该 算法有算法有p 阶精度。阶精度。第13页 欧拉法局部截断误差：欧拉法局部截断误差：Ri 主项主项/*leading term*/欧拉法含有欧拉法含有 1 阶精阶精度。度。第14页例例1:1:用欧拉公式求解初值问题用欧拉公式求解初值问题取步长取步长 。解解:应用应用EulerEuler公式于题给初值问题详细形式为：公式于题给初值问题详细形式为：其中其中 。计算结果列于下表：计算结果列于下表：第15页 第16页可用来检验近似解准确程度。可用来检验近似解准确程度。进行计算，数值解已到达了一定精度。进行计算，数值解已到达了一定精度。这个初值问题准确解为这个初值问题准确解为 ，从上表最终一列，我们看到取步长从上表最终一列，我们看到取步长第17页 欧拉公式改进：欧拉公式改进：隐式欧拉法隐式欧拉法/*implicit Euler method*/向后差商近似导数向后差商近似导数x0 x1)(,()(1101xyxfhyxy+第18页因为未知数因为未知数 yi+1 同时出现在等式两边，不能直接得同时出现在等式两边，不能直接得到，故称为到，故称为隐式隐式/*implicit*/欧拉公式，而前者称欧拉公式，而前者称为为显式显式/*explicit*/欧拉公式。欧拉公式。第19页普通先用显式计算一个初值，再普通先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法局部截断误差：欧拉法局部截断误差：即隐式欧拉公式含有即隐式欧拉公式含有 1 阶精度。阶精度。第20页 梯形公式梯形公式 /*trapezoid formula*/*trapezoid formula*/显、隐式两种算法显、隐式两种算法平均平均注：确有局部截断误差 ，即梯形公式即梯形公式含有含有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相同。迭代法，其迭代收敛性与欧拉公式相同。第21页 中点欧拉公式中点欧拉公式/*midpoint formula*/中心差商近似导数中心差商近似导数x0 x2x1假设假设 ,则能够导出则能够导出即中点公式含有即中点公式含有 2 阶精度。阶精度。第22页方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提升精度提升计算量大计算量大精度提升精度提升,显式显式多一个初值多一个初值,可能影响精度可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages?Do you think it possible?Well,call me greedy OK,lets make it possible.第23页 改进欧拉法改进欧拉法/*modified Eulers method*/Step 1:先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2:再将再将 代入代入隐式隐式梯形公式右边作梯形公式右边作校正校正，得到，得到1+iy第24页注注:此法亦称为此法亦称为预测预测-校正法校正法 /*predictor-corrector method*/*predictor-corrector method*/能够证实该算法能够证实该算法含有含有 2 阶精度阶精度，同时能够看到它，同时能够看到它是个是个单步单步递推格式，比隐式公式迭代求解过程递推格式，比隐式公式迭代求解过程简单简单。后面将看到，它。后面将看到，它稳定性高于稳定性高于显式欧拉法。显式欧拉法。第25页在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 应用改进欧拉法应用改进欧拉法,假如序列假如序列 收敛收敛,它极限便满足方程它极限便满足方程第26页所以，改进欧拉法公式含有所以，改进欧拉法公式含有 2 2 阶精度阶精度第27页例例2:2:用改进用改进Euler公式求解例公式求解例1中初值问题，中初值问题，取步长取步长 。解：解：对此初值问题采取改进对此初值问题采取改进EulerEuler公式，公式，其详细形式为其详细形式为计算结果列于下表：计算结果列于下表：第28页改进改进Euler法法Euler法法第29页经过计算结果得比较能够看出，改进经过计算结果得比较能够看出，改进Euler方法方法计算精度比计算精度比Euler方法要高。方法要高。第30页3 3 龙格龙格 -库塔法库塔法 /*/*Runge-KuttaRunge-Kutta Method*/Method*/建立高精度单步递推格式。建立高精度单步递推格式。单步递推法单步递推法基本思想基本思想是从是从(xi,yi)点出发，点出发，欧拉法及其各种变形所能到达最高精度为欧拉法及其各种变形所能到达最高精度为2 2阶。阶。以以某一斜率某一斜率沿直线到达沿直线到达 点。点。第31页 考查改进欧拉法，能够将其改写为：考查改进欧拉法，能够将其改写为：斜率斜率一定取一定取K1 K2 平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？第32页将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+首先希望能确定系数首先希望能确定系数 1、2、p，使得到算法格式，使得到算法格式有有2阶精度，即在阶精度，即在 前提假设下，使得前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开第33页Step 2:将将 K2 代入第代入第1式，得到式，得到Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点点泰勒泰勒展开作比较展开作比较第34页要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。全部满足上式格式统称为。全部满足上式格式统称为2阶龙阶龙格格-库塔格式库塔格式。注意到，注意到，就是改进欧拉法。就是改进欧拉法。第35页 为取得更高精度，应该怎样深入推广？为取得更高精度，应该怎样深入推广？).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待定系数，确定这均为待定系数，确定这些系数步骤与前面相同。些系数步骤与前面相同。第36页v考虑一阶常微分方程初值问题考虑一阶常微分方程初值问题第37页将区域将区域a,ba,b进行分划：进行分划：第38页 若若第39页则则第40页n n级显式级显式Runge-KuttaRunge-Kutta方法方法第41页二阶二阶Runge-KuttaRunge-Kutta方法方法取取n=2n=2记记第42页由此得由此得第43页v其次为使局部截断误差为为使局部截断误差为 ,应取应取第44页改进改进EulerEuler方法方法v取取第45页中点方法中点方法取取 第46页二阶二阶HeunHeun方法方法取取 第47页二级二级Runge-KuttaRunge-Kutta方法不超出二阶方法不超出二阶v记记则则第48页所以局部截断误差只能到达所以局部截断误差只能到达第49页三级三级Runge-KuttaRunge-Kutta方法方法取取n=3n=3第50页记记第51页第52页v又因为又因为第53页所以要使局部截断误差为所以要使局部截断误差为 ，必须，必须 第54页KuttaKutta方法方法取取第55页三阶三阶HeunHeun方法方法v取取第56页三级三级Runge-KuttaRunge-Kutta方法不超出三阶方法不超出三阶v完全类似于二级完全类似于二级Runge-KuttaRunge-Kutta方法分析方法分析只能到达只能到达三级三级Runge-KuttaRunge-Kutta方法局部截断误差方法局部截断误差将将 和和 都展开到都展开到 项项,易证易证第57页 四级四级R-KR-K方法方法取取n=4n=4第58页第59页局部截断误差为O(h5)四阶四阶经典龙格经典龙格-库塔法库塔法 /*Classical Runge-Kutta Method*/：第60页附注：附注：二阶二阶Runge-KuttaRunge-Kutta方法局部截断误差方法局部截断误差 只能到达只能到达 五阶五阶Runge-KuttaRunge-Kutta方法局部截断误差方法局部截断误差 只能到达只能到达 四阶四阶Runge-KuttaRunge-Kutta方法局部截断误差方法局部截断误差 只能到达只能到达 三阶三阶Runge-KuttaRunge-Kutta方法局部截断误差方法局部截断误差 只能到达只能到达 第61页注：注：龙格龙格-库塔法主要运算在于计算库塔法主要运算在于计算 值，即计算值，即计算 值。值。Butcher 于于1965年给出了计算量与可到达最高年给出了计算量与可到达最高精度阶数关系：精度阶数关系：753可到达最高精度可到达最高精度642每步须算每步须算Ki 个数个数 因为龙格因为龙格-库塔法导出基于泰勒展开，故精库塔法导出基于泰勒展开，故精太好解，最好采取太好解，最好采取低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数光滑性影响。对于光滑性不度主要受解函数光滑性影响。对于光滑性不第62页4 4 单步方法收敛性与稳定性单步方法收敛性与稳定性 /*/*Convergency and StabilityConvergency and Stability*/*/收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定 x=xi=x0+i h，当当 h0(同同时时 i )时时有有 yi y(xi)，则则称称该该算算法法是是收敛收敛。例：例：就初值问题就初值问题 考查欧拉显式格式收考查欧拉显式格式收敛性。敛性。第63页解：解：该问题准确解为该问题准确解为 欧拉公式为欧拉公式为对任意固定对任意固定 x=xi=i h，有，有 第64页 稳定性稳定性/*Stability*/例：例：考查初值问题考查初值问题 在区间在区间0,0.5上解上解.分别用欧拉显、隐式格式和改进欧拉格式分别用欧拉显、隐式格式和改进欧拉格式计算数值解。计算数值解。0.00.10.20.30.40.5准确解准确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.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?!An Engineer complains:Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!第65页若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生误误差差在在以以后后计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定/*absolutely stable*/。普通分析时为简单起见，只考虑普通分析时为简单起见，只考虑试验方程试验方程/*test equation*/常数，能够常数，能够是复数是复数第66页我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 绝对稳定绝对稳定区域比区域比 B 大大。当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在时，将某算法应用于上式，并假设只在初值产生误差初值产生误差 ，则若此误差以后逐步衰减，则若此误差以后逐步衰减,就称该算法相对于就称该算法相对于 绝对稳定，绝对稳定，全体组成全体组成h h=h绝对稳定区域绝对稳定区域。第67页例：例：考查显式欧拉法考查显式欧拉法0-1-2ReImg由此可见，要确保初始误差由此可见，要确保初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：第68页例：例：考查隐式欧拉法考查隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：普通来说，隐式欧拉法绝对稳定性比同阶显普通来说，隐式欧拉法绝对稳定性比同阶显式法好。式法好。第69页例：例：隐式龙格隐式龙格-库塔法库塔法其中其中2阶方法阶方法 绝对稳定绝对稳定区域为区域为0ReImg无条件稳定无条件稳定第70页而而显式显式 1 4 阶方法绝对稳定区域为阶方法绝对稳定区域为k=1k=2k=3k=4-1-2-3-123ReImg第71页