《最新常微分方程数值解aPPT课件.ppt》由会员分享，可在线阅读，更多相关《最新常微分方程数值解aPPT课件.ppt（24页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、常微分方程数值解常微分方程数值解a2 欧拉方法欧拉方法/*Eulers Method*/欧拉公式：欧拉公式：x0 x1向前差商近似导数向前差商近似导数记为记为在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的的截截断断误误差差 Ri=y(xi+1)yi+1 称称为为局局部部截截断断误误差差/*local truncation error*/。若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。欧拉法的局部截断误差：欧拉法的局部截断误差：欧拉法具有欧拉法具有 1 阶精度。阶精度。其中其中 i
2、(i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算
3、f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。4 收敛性与稳定性收敛性与稳定性 /*Convergency and Stability*/收敛性收敛性/*Convergency*/若若某某算算法
4、法对对于于任任意意固固定定的的 x=xi=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 稳定性稳定性/*Stability*/例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确
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 7 An Engineer complains:Math theorems ar
6、e so unstable that a small perturbation on the conditions will cause a crash on the conclusions!若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定的的/*absolutely stable*/。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程/*test equation*/常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在
7、初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该算法，则若此误差以后逐步衰减，就称该算法相对于相对于 绝对稳定绝对稳定，的全体构成的全体构成绝对稳定区域绝对稳定区域。我们。我们称称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。h h=h例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210注：注：一般来说，隐式欧拉法的
8、绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。5 线性多步法线性多步法 /*Multistep Method*/用用若干若干节点处的节点处的 y 及及 y 值的值的线性组合线性组合来近似来近似y(xi+1)。).(.110111101kikiiikikiiiffffhyyyy +=其通式可写为：其通式可写为：当当 1 0 时，为时，为隐式公式隐式公式;1=0 则为则为显式公式显式公式。基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通过，则可通过 近似近似y(xi+1)。而。
9、而选用不同近似式选用不同近似式 Ik，可得到不同的计算公式，可得到不同的计算公式。亚当姆斯显式公式亚当姆斯显式公式/*Adams explicit formulae*/利用利用k+1 个节点上的被积函数值个节点上的被积函数值 构造构造 k 阶牛顿阶牛顿后后插插多项式多项式 ,有有Newton插值余项插值余项/*显式计算公式显式计算公式*/局部截断误差为：局部截断误差为：例：例：k=1 时有时有注：注：一般有一般有 ，其中，其中Bk 与与yi+1 计算公式中计算公式中 fi,fi k 各项的各项的系数系数均可查表得到均可查表得到。10123kfifi 1fi 2fi 3Bk常用的是常用的是 k=
10、3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式 亚当姆斯隐式公式亚当姆斯隐式公式/*Adams implicit formulae*/利用利用k+1 个节点上的被积函数值个节点上的被积函数值 fi+1,fi,fi k+1 构造构造 k 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项式完全类似地可得到一系列隐式公式隐式公式，并有，并有 ，其中，其中 与与 fi+1,fi,fi k+1 的系数亦可查表得到。的系数亦可查表得到。10123kfi+1fifi 1fi 2Bk常用的是常用的是 k=3 的的4阶亚当姆斯隐式公式阶亚当姆斯隐式公式小于小于Bk较同阶显较同阶显
11、式式稳定稳定 亚当姆斯预测亚当姆斯预测-校正系统校正系统 /*Adams predictor-corrector system*/Step 1:用用Runge-Kutta 法法计算前计算前 k 个初值；个初值；Step 2:用用Adams 显式显式计算计算预测预测值；值；Step 3:用同阶用同阶Adams 隐式隐式计算计算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必须相同。式的精度必须相同。通常用通常用经典经典Runge-Kutta 法法配合配合4阶阶Adams 公式公式。Hey!Look at the local truncation error of the explici
12、t and implicit Adams methods:and Dont you think theres something you can do?4阶阶Adams隐式公式的截断误差为隐式公式的截断误差为4阶阶Adams显式公式的截断误差为显式公式的截断误差为当当 h 充分小充分小时，可近似认为时，可近似认为 i i，则：，则：Predicted value pi+1Modified value mi+1Corrected value ci+1Modified final value yi+1外推技术外推技术/*extrapolation*/基于泰勒展开的构造法基于泰勒展开的构造法).(.
13、110111101kikiiikikiiiffffhyyyy +=将通式中的右端各项将通式中的右端各项 yi 1,yi k;fi+1,fi 1,fi k 分别在分别在 xi 点作点作泰勒展开泰勒展开，与精确解，与精确解 y(xi+1)在在 xi 点的泰勒展开作点的泰勒展开作比较比较。通过令。通过令同同类项系数相等类项系数相等，得到足以确定待定系数，得到足以确定待定系数 0,k;1,0,k 的等式，则可构造出线性多步法的的等式，则可构造出线性多步法的公式。公式。例：例：设设)(3322110221101 +=iiiiiiiiyyyyhyyyy 确定式中待定系数确定式中待定系数 0,1,2,0,1
14、,2,3,使得公式具有使得公式具有4阶阶精度。精度。解：解：/*y(xi)=yi*/个未知数个未知数个方程个方程75 令令 1=2=0Adams 显式显式公式公式 以以 y i+1 取代取代 y i 1，并取，并取 1=2=0Adams 隐式隐式公式公式 以以 yi 3 取代取代 y i 3，则可导出另一组，则可导出另一组4 阶显式算法，其阶显式算法，其中包含了著名的中包含了著名的米尔尼米尔尼/*Milne*/公式公式其局部截断误差为其局部截断误差为注：注：上式也可通过上式也可通过数值积分数值积分导出，即将导出，即将 在区间在区间 上积分，得到上积分，得到 再再过过 做做 f 的插值多项式即可
15、。的插值多项式即可。取取 1=1,2=0得到得到辛甫生辛甫生/*Simpson*/公式公式与与Milne 公式匹配使用公式匹配使用 辛甫生辛甫生/*Simpson*/公式公式在区间在区间xi 1,xi+1上积分，并用上积分，并用Simpson数值积分数值积分公式来近似积公式来近似积分项，亦可得此分项，亦可得此Simpson公式。公式。Milne-Simpson 系统的缺点是系统的缺点是稳定性差稳定性差，为改善稳定性，为改善稳定性，考虑另一种隐式校正公式：考虑另一种隐式校正公式：要求公式具有要求公式具有4 阶精度。通过泰勒展开，可得到阶精度。通过泰勒展开，可得到 个等式，个等式，从中解出从中解出 个未知数，则有个未知数，则有 个自由度。个自由度。561取取 1=1 得得Simpson 公式公式哈明哈明/*Hamming*/用用 1 的不同数值进行试验，发现当的不同数值进行试验，发现当 1=0 时，公式的稳定性较好，即：时，公式的稳定性较好，即：其局部截断误差为其局部截断误差为注：注：哈明公式不能用数值积分方法推导出来。哈明公式不能用数值积分方法推导出来。