《数值计算》课程教学课件（讲稿）第8章 常微分方程的数值解法（2/2）

m教材《数值分析》陈晓江主编武汉理工大学出版社Y参考书目《数值分析》李庆扬、王能超、易大义编清华大学出版社施普林格出版社《科学和工程计算基础》施妙根、顾丽珍编清华大学出版社上页下页返圆

上页 下页 返回 教材 《数值分析》 陈晓江主编 武汉理工大学出版社  参考书目 《数值分析》 李庆扬、王能超、易大义编 清华大学出版社 施普林格出版社 《科学和工程计算基础》 施妙根、顾丽珍编 清华大学出版社

第八章常微分方程的数值解法第一节问题的提出第二节欧拉法第三节龙格-库塔法第四节单步法的收敛性与稳定性第五节线性多步法第六节一阶方程组和高阶方程第七节边值问题的数值解法上页下页返圆

上页 下页 返回 第八章 常微分方程的数值解法 第一节 问题的提出 第二节 欧拉法 第三节 龙格-库塔法 第四节 单步法的收敛性与稳定性 第五节 线性多步法 第六节 一阶方程组和高阶方程 第七节 边值问题的数值解法

$4单步法的收敛性与稳定性一.收敛性定义若某算法对于任意固定的x=x,=xo+ih，当h>0(同时i→o)时有yi→(x;)，则称该算法是收敛的。y'=2y考察欧拉显式格式的收敛性。例：就初值问题y(0) = yo解： 该问题的精确解为 (x)= yoeax欧拉公式为yi+i=y;+hay;=(1+ah)ylim(1+ ah)/ah=e+h->0对任意固定的x=x;=ih，有J; = o(1 + 2h)*; /h上页= yo[(1+ 2h)/ahjai → oeai = y(x,)下页返圆

上页 下页 返回 一. 收敛性 定义 若某算法对于任意固定的 x = xi = x0 + i h，当 h0 ( 同时 i  ) 时有 yi  y( xi )，则称该算法是收敛的。 例：就初值问题 考察欧拉显式格式的收敛性。        0 y(0) y y y 解：该问题的精确解为 x y x y e  0 ( )  欧拉公式为 i i i i y y h y (1 h) y 1       0 y (1 h) y i i    对任意固定的x = xi = i h ，有 x h i i y y h / 0  (1   ) h e h h      1/ 0 lim(1 ) xi y e   0  §4 单步法的收敛性与稳定性 ( ) xi  y h xi y h    [(1   ) ] 1/ 0

二.稳定性y(x)=-30y(x) ,在区间[0,0.5]上的解。例：考察初值问题y(0)=1分别用欧拉显、隐式格式和改进的欧拉格式计算数值解J=e-30x节点X欧拉显式欧拉隐式精确解改进欧拉法0.01.00001.00001.00001.00000.12.5000x10-14.9787×10-2-2.00002.50000.26.2500x10-22.4788x10-34.00006.25000.31.5625x10-21.5626x1011.2341x10-4-8.00000.43.9063x10-33.9063x1016.1442x10-61.6000x1010.53.0590x10-79.7656x10-49.7656x101-3.2000×101上页下页返圆

上页 下页 返回 二. 稳定性 例：考察初值问题 在区间[0, 0.5]上的解。 分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。        (0) 1 ( ) 30 ( ) y y x y x 0.0 0.1 0.2 0.3 0.4 0.5 节点 xi 欧拉显式 欧拉隐式 改进欧拉法 精确解 x y e 30  1.0000 2.0000 4.0000 8.0000 1.6000101 3.2000101 1.0000 2.5000101 6.2500102 1.5625102 3.9063103 9.7656104 1.0000 2.5000 6.2500 1.5626101 3.9063101 9.7656101 1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107

定义若某算法在计算过的误差在以后的计入为复数。为保算中都逐步衰减，则证稳定性，通常取入的实部小于0一般分析时为简单起见，只广开王y'=ay当步长取为h时，将某算法应用于上式，并假设只在初值产生误差&。=Jy。-J。，则若此误差以后逐步衰减，就称该算法相对于h=ah绝对稳定，h的全体构成绝对稳定区域。我们称算法A比算法B稳定，就是指A的绝对稳定区域比B的大。上页下页返圆

上页 下页 返回 定义 若某算法在计算过程中任一步产生的误差在以后的计 算中都逐步衰减，则称该算法是绝对稳定的。 一般分析时为简单起见，只考虑试验方程 y   y λ为复数。为保 证稳定性，通常 取λ的实部小于0 当步长取为h 时，将某算法应用于上式，并假设只在初值 产生误差 ，则若此误差以后逐步衰减，就称该 算法相对于 绝对稳定， 的全体构成绝对稳定区域。 我们称算法A 比算法B 稳定，就是指A 的绝对稳定区域比B 的大。 0 0 0   y  y h   h h

例：考察显式欧拉法yi+r=y,+hay,=(l+h)iyoImgCo = yo- y。 Ji+i=(1+h)i+l joRe8i+1 = yi+1 - Ji+1 =(1+h)i+l o由此可见，要保证初始误差&以后逐步衰减，h=ah必须满足：|1+h|1注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式上页法的好。下页返圆

上页 下页 返回 例：考察显式欧拉法 0 1 1 y y h y (1 h) y i i i i        0 0 0   y  y 0 1 1 y (1 h) y i i     0 1 1 1 1  (1 )          i i yi yi h 由此可见，要保证初始误差0 以后逐步衰减， 必须满足： h   h |1 h|  1 -2 -1 0 Re Img 例：考察隐式欧拉法 i1  i  i1 y y h y i i y h y          1 1 1 0 1 1 1 1             i i h 可见绝对稳定区域为： |1h|  1 0 1 2 Re Img 注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式 法的好

Ji+i =y; +h[a,K, +... +amKm]例：隐式龙格-库塔法K, =f(x, +a,h, y, +β,hK, +... +βjmhKm)(j=l,...,m)Yi+I = y; +hK,hh其中2阶方法的绝对稳定区域为K,=f(x,+,y +"K.22个Img而显式1~4阶方法的绝对稳定区域为个Imgk=4330Rek=2k=1>-3-2-1Re无条件稳定上页下页返回

上页 下页 返回 例：隐式龙格-库塔法                 ( 1, . , ) ( , . ) [ . ] 1 1 1 1 1 j m K f x h y hK hK y y h K K j i j i j jm m i i m m      而显式 1~ 4 阶方法的绝对稳定 区域为           ) 2 , 2 ( 1 1 1 1 K h y h K f x y y hK i i i i 其中2阶方法 的绝对稳定区域为 0 Re Img k=1 k=2 k=3 k=4 -3 -2 -1 - - - 1 2 3 Re Img 无条件稳定

85线性多步法当β-时，为隐式公式;两步欧拉公式：Ji+1=y用若干节点处的y及βI=0则为显式公式。其通式可写为：Yi+ =αy, +αi yi- +.+αyi-k +h(β-. fi+ +βofi + β fi- +.+ β fi-k)>基于数值积分的构造法将'=f(x,y)在[x;,xitil 上积分，得到y(xi+1) - y(x,) = (x* f(x, y(x)dx只要近似地算出右边的积分Ik~f(x,y(x))dx，则可通过yi+1=y;+I近似y(xi+1)。而选用不同近似式Ik，可得到不上页同的计算公式。下页返圆

上页 下页 返回 用若干节点处的y 及 y’ 值的线性组合来近似y(xi+1)。 . ( . ) i 1 0 i 1 i 1 k i k 1 i 1 0 i 1 i 1 k i k y y y y h f f f f                     其通式可写为： ( , ) j j j f  f x y  基于数值积分的构造法 将 y  f (x, y) 在 [xi , xi1 ] 上积分，得到      1 ( ) ( ) ( , ( )) 1 i i x x y xi y xi f x y x d x 只要近似地算出右边的积分 ，则可通 过 近似y(xi+1)。而选用不同近似式Ik，可得到不 同的计算公式。    1 ( , ( )) i i x x I k f x y x dx i i k y  y  I 1 §5 线性多步法 两步欧拉公式： yi1  yi1  2h 当 f (xi , yi ) i  1, . ,n 1 10 时，为隐式公式; 1=0 则为显式公式

亚当姆斯显式公式利用k+1个节点上的被积函数值fi,fi-,,fi-k构造k阶牛顿后插多项式N,(x,+th)，te[0,1],有[*" f(x,y(x)dx=["N,(x, +th)hdt+'R,(x, +th)hdtYi+ = y;+hf"'N,(x, +th)dt牛顿插值余项局部截断误差为：R;=(xi+1)-i+1=h[,R(x,+th例: k-1 时有 N;(x; +th)=f;+t(f;-fi-1)hYi1 = y +hf,lf, +t(f; - fi-)Idt = y; +=(3f - fil)2R, =hf"d"f(s,y(5) 15上页h'y"()th(t +1)h dt =dx?02!12下页返圆

上页 下页 返回  亚当姆斯显式公式 利用k+1 个节点上的被积函数值 构造 k 阶牛顿 后插多项式 , 有 i i i k f f f   , , . , 1 N (x  th) , t [0, 1] k i         1 0 1 0 ( , ( ) ) ( ) ( ) 1 f x y x d x N x th hd t R x th hd t k i x x k i i i 牛顿 插值余项 y y h N x th dt i i k i ( ) 1 0 1      局部截断误差为：        1 0 1 1 R y(x ) y h R (x th)d t i i i k i 例：k=1 时有 ( ) ( ) 1 i   i  i  i1 N x th f t f f            1 0 1 1 1 (3 ) 2 [ ( )] i i i i i i i i f f h y y h f t f f dt y th t h d t d x d f y R h x x i ( 1) 2! 1 ( , ( )) 1 0 2 2      ( ) 12 5 3 i  h y 

注：一般有R, =B,hk+2 y(k+2)(5)，其中B,与yi+i计算公式中fi....fi-各项的系数均可查表得到。kBkJifi-1 fi-2 Ji-31-252013-21238231652121212552515937324242472024常用的是k=3的4阶亚当姆斯显式公式上页h(55f -59 f- +37f-, -9 f-3)Ji+1 =Y, +下页24返圆

上页 下页 返回 注：一般有 ，其中Bk 与yi+1 计算公式 中 fi , ., fik 各项的系数均可查表得到。 ( ) 2 ( 2) i k k i k R B h y     0 1 1 2 3 k 2 1 2 3 12 23 24 55 2 1  12 16  24 59  12 5 24 37 24 9  12 5 8 3 720 251 fi fi1 fi2 fi3 . Bk . . . . . . . 常用的是 k = 3 的4阶亚当姆斯显式公式 (55 59 37 9 ) 24 i1  i  i  i1  i2  i3 f f f f h y y