数值分析(研究生)第七章常微分方程的数值解法二

1、第七章第七章 常微分方程的数值解法（二）常微分方程的数值解法（二）第四节第四节 线性多步法线性多步法第五节第五节 单步法的收敛性与稳定性单步法的收敛性与稳定性第六节第六节 一阶方程组和高阶方程一阶方程组和高阶方程第七节第七节 边值问题的数值解法边值问题的数值解法4 线性多步法线性多步法用用若干若干节点处的节点处的 y 及及 y 值的值的线性组合线性组合来近似来近似y(xi+1).).(.110111101kikiiikikiiiffffhyyyy- - -+ +- - - -+ + + + + + + + + += =b bb bb bb ba aa aa a其通式可写为：其通式可写为：),(

2、jjjyxff = =当当 b b- -1 0 时，为时，为隐式公隐式公式式; b b- -1=0 则为则为显式公式显式公式. 基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxfy =,1+iixx + += =- -+ +1)(,()()(1iixxiidxxyxfxyxy只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通，则可通过过 近似近似y(xi+1) .而而选用不同近似式选用不同近似式 Ik，可得到不同，可得到不同的计算公式的计算公式. + + 1)(,(iixxkdxxyxfIkiiIyy+ += =+ +1 亚当姆斯显式公式亚当姆斯显

3、式公式利用利用k+1 个节点上的被积函数值个节点上的被积函数值 构造构造 k 阶牛顿阶牛顿后插后插多项式多项式 , 有有kiiifff- - -,.,11, 0, )( + +thtxNik + + + += =+ +1010)()()(,(1dthhtxRdthhtxNdxxyxfikxxikiiNewton插值余项插值余项dthtxNhyyikii)(101+ + += = + +/* 显式计算公式显式计算公式 */局部截断误差为：局部截断误差为： + += =- -= =+ + +1011)()(dthtxRhyxyRikiii例例1 k=1 时有时有)()(11- - -+ += =

4、+ += =+ +iiiiiifftfftfhtxN - - -+ +- -+ += =- -+ + += =10111)3(2)(iiiiiiiiffhydtfftfhyydththtdxyfdhRxxi)1(!21)(,(1022+ += = )(1253iyh = =注：注：一般有一般有 ，其中，其中Bk 与与yi+1 计算公式计算公式中中 fi , , fi- -k 各项的各项的系数系数均可查表得到均可查表得到 . )()2(2ikkkiyhBR + + += =10123k21231223245521-1216-2459-1252437249-12583720251fifi- -1f

5、i- -2fi- -3BkMisprint on p.106常用的是常用的是 k = 3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式)9375955(243211- - - -+ +- -+ +- -+ += =iiiiiiffffhyy 亚当姆斯隐式公式亚当姆斯隐式公式利用利用k+1 个节点上的被积函数值个节点上的被积函数值 fi+1 , fi , , fi- -k+1 构造构造 k 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项式完全类似地可得到一系列隐式公式隐式公式，并有，并有 ，其中，其中 与与 fi+1 , fi , , fi- -k+1 的系数亦可

6、查表得到的系数亦可查表得到.)()2(2ikkkiyhBR + + += =kB10123k21- -21125249211282419121- -245- -241121- -241- -72019- -fi+1fifi- -1fi- -2Bk常用的是常用的是 k = 3 的的4阶亚当姆斯隐式公式阶亚当姆斯隐式公式)5199(242111- - -+ + + +- -+ + += =iiiiiiffffhyy小于小于Bk较同阶显较同阶显式式稳定稳定 亚当姆斯预测亚当姆斯预测-校正系统校正系统 Step 1: 用用Runge-Kutta 法法计算前计算前 k 个初值；个初值；Step 2: 用

7、用Adams 显式显式计算计算预测预测值；值；Step 3: 用同阶用同阶Adams 隐式隐式计算计算校正校正值值.注意：注意：三步所用公三步所用公式的精度必须相同。式的精度必须相同。通常用通常用经典经典Runge-Kutta 法法配合配合4阶阶Adams 公式公式.4阶阶Adams隐式公式的截断误差为隐式公式的截断误差为)(72019)()5(511iiiyhyxy - -= =- -+ + +4阶阶Adams显式公式的截断误差为显式公式的截断误差为)(720251)()5(511iiiyhyxy = =- -+ + +当当 h 充分小充分小时，可近似认为时，可近似认为 i i ，则：，则：

8、19251)()(1111- - - - -+ + + + +iiiiyxyyxy)(270251)(1111+ + + + +- -+ + iiiiyyyxy)(27019)(1111+ + + + +- - - iiiiyyyxyPredicted value pi+1Modified value mi+1Corrected value ci+1Modified final value yi+1外推技术外推技术 /* extrapolation */.).,(),(2211尽可能高的阶尽可能高的阶，使该公式具有，使该公式具有、展开原理确定展开原理确定试用试用其中其中给出一个三步显式格式给出

9、一个三步显式格式例例b ba ab ba aTayloryxffffhyynnnnnnn= =+ + += =- - -+ +5 收敛性和稳定性收敛性和稳定性一、一、收敛性收敛性定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的的. 例例3 就初值问题就初值问题 考察欧拉显式格式的收敛性考察欧拉显式格式的收敛性. = = = 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)(= =尤拉公式为尤拉公式为iiiiyhyhyy)

10、1 (1 + += =+ += =+ +0)1 (yhyii + += =对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0+ += =+ += =ehhh=+ /10)1(lim)(0ixxyeyi= = .1)0(03*3的的梯梯形形法法收收敛敛证证明明：解解初初值值问问题题例例 = = =+ + yyy0.00.10.20.30.40.5精确解精确解改进尤拉法改进尤拉法 尤拉隐式尤拉隐式尤拉显式尤拉显式 节点节点 xixey30-=二、二、 稳定性稳定性例例4 考察初值问题考察初值问题 在区间在区间0, 0.5上的解上的解.分别用

11、欧拉显、隐式格式和改进的欧拉格式计算数值解分别用欧拉显、隐式格式和改进的欧拉格式计算数值解. = =- -= = 1)0()(30)(yxyxy 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.144

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

13、稳定区域.我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大.000yy - -= = h h= =h例例5 考察显式欧拉法考察显式欧拉法011)1(yhyhyyiiii+ + + += =+ += = 000yy - -= = 011)1(yhyii+ + + += =01111)1( + + + + + += =- -= =iiiihyy由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh = =1|1| + + h0-1-2ReImg例例6 考察隐式欧拉法考察隐式欧拉法

14、11+ + + += =iiiyhyy iiyhy - -= =+ +11101111 + + + - -= =iih可见绝对稳定区域为：可见绝对稳定区域为：1|1| - -h210ReImg注：注：一般来说，隐式尤拉法的绝对稳定性比同阶的显式一般来说，隐式尤拉法的绝对稳定性比同阶的显式法的好法的好.例例7 隐式龙格隐式龙格-库塔法库塔法 = =+ + + + += =+ + + += =+ +),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmiib bb ba a 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 + + += =+ +=

15、 =+ +)2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定6 微分方程组与高阶方程微分方程组与高阶方程一、一、一阶微分方程组一阶微分方程组IVP的一般形式为：的一般形式为：= = = = )(,.),(,()(.)(,.),(,()(1111xyxyxfxyxyxyxfxymmmm初值初值0002020101)(,.,)(,)(mmyxyyxyyxy= = = =将问题记作将问题记作向量形式向量形式，令：，令： = = = = = =001011.,.,

16、.mmmyyyfffyyy = = = 00)(),()(yxyyxfxy前述所有公式皆前述所有公式皆适用于向量形式适用于向量形式.二、二、高阶微分方程高阶微分方程= = = = = = =- - - -10)1(1000)1()()(,.,)(,)(),.,(nnnnaxyaxyaxyyyyxfy化作化作一阶微分方程组一阶微分方程组求解求解.引入新变量引入新变量)1(21,.,- -= = = = =nnyyyyyy = = = = = = - -),.,(.1121nnnnyyxfyyyyy初值条件为：初值条件为：10102001)(.)()(- -= = = =nnaxyaxyaxy).

17、2 . 0()2 . 0(2 . 0)(2)0(, 1)0(1)(8yyxxyyyyxyyyxyKuttaRunge = = = = = = = = =+ + + +- - - -、处处的的函函数数值值及及其其导导数数值值在在的的解解初初值值问问题题格格式式求求二二阶阶常常微微分分方方程程用用四四阶阶例例7 边值问题的数值解法边值问题的数值解法2 阶常微分方程边值问题阶常微分方程边值问题 = = = = = b ba a)(,)(),(),(byaybaxyyxfy 打靶法打靶法先猜测一个初始斜率先猜测一个初始斜率 y (a) = s，通过解初值，通过解初值问题问题 = = = = = = sayaayyyxfy)()(),(y(b) = (s)找出找出s*使得使得 (s*) = b b，即把问，即把问题转化为求方程题转化为求方程 (s) - - b b = 0 的根的根.yx0aby x( )b b斜率斜率 = s0 ()s0斜率斜率 = s1 ()s1每计算一个每计算一个 (s) 都必须解一个都必须解一个ODE. 有限差分法有限差分法将求解区间将求解区间a, b 等分为等分为N 份，取节点份，取