库塔格式注意到就是改进的欧拉法课件

第六章第六章 常微分方程数值解常微分方程数值解/*Numerical Methods for Ordinary Differential Equations*/考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题/*Initial-Value Problem*/:只要只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与，即存在与 x,y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a,b 上的上的 y1(x)和和 y2(x)都成立，则上述都成立，则上述IVP存在唯一解存在唯一解。要计算出解函数要计算出解函数 y(x)在一系列节点在一系列节点 a=x0 x1 xn=b 处的近似值处的近似值节点间距节点间距 为步长，通常采用为步长，通常采用等距节等距节点点，即取，即取 hi=h(常数常数)。第六章 常微分方程数值解 考虑一阶常微分方程的初值问题11 欧拉方法欧拉方法/*Eulers Method*/欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数记为记为在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的的截截断断误误差差 Ri=y(xi+1)yi+1 称称为为局局部部截截断断误误差差/*local truncation error*/。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。欧拉法的局部截断误差：欧拉法的局部截断误差：欧拉法具有欧拉法具有 1 阶精度。阶精度。Ri 的的主项主项/*leading term*/亦称为亦称为欧拉折线法欧拉折线法/*Eulers polygonal arc method*/1 欧拉方法/*Eulers Method*/2 欧拉公式的改进：欧拉公式的改进：隐式欧拉法隐式欧拉法/*implicit Euler method*/向后差商近似导数向后差商近似导数x0 x1)(,()(1101xyxfhyxy+)1,.,0(),(111=+=+niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显式显式/*explicit*/欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。Hey!Isnt the leading term of the local truncation error of Eulers method?Seems that we can make a good use of it 欧拉公式的改进：隐式欧拉法/*implicit3 梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法的显、隐式两种算法的平均平均注：注：的确有局部截断误差的确有局部截断误差 ，即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。中点欧拉公式中点欧拉公式/*midpoint formula*/中心差商近似导数中心差商近似导数x0 x2x1假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法/*double-step method*/，而前面的三种算法都是，而前面的三种算法都是单步法单步法/*single-step method*/。梯形公式/*trapezoid formula*/4方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高,显式显式多一个初值多一个初值,可能影响精度可能影响精度 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.方 法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性5 改进欧拉法改进欧拉法/*modified Eulers method*/Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法/*predictor-corrector method*/。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。改进欧拉法/*modified Eulers me62 龙格龙格-库塔法库塔法/*Runge-Kutta Method*/建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一某一斜率斜率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形点。欧拉法及其各种变形所能达到的最高精度为所能达到的最高精度为2阶阶。考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？2 龙格-库塔法/*Runge-Kutta Me7首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入第代入第1式，得到式，得到首先希望能确定系数 1、2、p，使得到的算法格式有2阶精8Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点的点的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？Step 3:将 yi+1 与 y(xi+1)在 x9其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：其中i (i=1,m)，i (i 10注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。注：753可达到的最高精度642每步须算Ki 的个数 由于113 收敛性与稳定性收敛性与稳定性 /*Convergency and Stability*/收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定的的 x=xi=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 3 收敛性与稳定性/*Convergency and12 稳定性稳定性/*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?!An Engineer complains:Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!稳定性/*Stability*/例：考察初值问题 13若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定的的/*absolutely stable*/。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程/*test equation*/常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定，的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。h h=h定义若某算法在计算过程中任一步产生的误差在以后的计算中14例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。例：考察显式欧拉法由此可见，要保证初始误差0 以后逐步衰减15例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定例：隐式龙格-库塔法而显式 1 4 阶方法的绝对稳定区域为16期末考试：期末考试：2.第三章第三章 非线性方程求根非线性方程求根：二分法、迭代法、牛顿法和弦截法：二分法、迭代法、牛顿法和弦截法要求：根的存在，公式，收敛性条件的判别要求：根的存在，公式，收敛性条件的判别1.第二章第二章 解线性方程组的直接法：掌握解线性方程组的直接法：掌握Gauss消元法进行到底消元法进行到底 的条件，矩阵三角分解定理的条件和结论，向量和的条件，矩阵三角分解定理的条件和结论，向量和 矩阵的范数，方程组的条件数与病态方程组的求解矩阵的范数，方程组的条件数与病态方程组的求解解线性方程组的迭代法：雅可比迭代法，高斯解线性方程组的迭代法：雅可比迭代法，高斯-赛德尔赛德尔 迭代法；要求：求解公式，收敛条件。迭代法；要求：求解公式，收敛条件。期末考试：2.第三章 非线性方程求根：二分法、迭代法、牛174.函数插值：拉格朗日插值，牛顿插值，埃米尔特插值函数插值：拉格朗日插值，牛顿插值，埃米尔特插值 要求：插值公式，余项公式要求：插值公式，余项公式 5.数值积分：插值型求积公式（矩形、梯形、辛普生、数值积分：插值型求积公式（矩形、梯形、辛普生、龙贝格）。公式和误差，代数精确度的概念。龙贝格）。公式和误差，代数精确度的概念。6.常微分方程数值解：单步法（欧拉、改进欧拉，龙格常微分方程数值解：单步法（欧拉、改进欧拉，龙格-库塔）库塔）差分法。要求：求解公式和误差阶。差分法。要求：求解公式和误差阶。4.函数插值：拉格朗日插值，牛顿插值，埃米尔特插值5.数18