数值分析Ch9常微分方程数值解法

数 值 分 析 Computational Method Chapter 9 常微分方程数值解法 9.常微分方程数值解法 9．1 引言 • 一阶常微分方程初值问题可为:      00 )( ),( yxy yxfy 若函数 适当光滑，且关于y满足李普希兹 (Lipschitz)条件 ),( yxf 2121 ),(),( yyLyxfyxf  则上述问题存在唯一连续可微的解 y(x)。 数值方法基本思想：在初值问题的求解区间[a,b] 上插入若干个节点： 1 iii xxh 称为步长 数值解法的特点:采用步进式或递推式算法. 由初始值 y(x0)推导出x1的近似解 y1，再由 y1计算出 y2，y3，，直到求出 yn为止。沿节点次序一步一步 地向前推进. 总假定等步长 hhi  niihxxi ,,2,1,00  数值解法: nxxx ,,, 10  将常微分方程这一连续模型离散化,求 bxxxxxa nn  1210  处的近似值 .,,, 10 nyyy  得初值问题的解函数y(x)在离散节点 由yi-k,yi-k+1, ,yi-2 , yi-1得到yi 称为多步法。或k步法 仅由yi-1得到yi 称为单步法。 所求出的 yi (i=1,2,,n)称为初值问题的数值解。 9.2.简单的数值方法与基本概念 • 1.微分法: h xyxy xy nnn )()( )( 1    ),(),()( 1 nnnnnnn yxhfyyyxfxy   欧拉(显式) h xyxy xy nnn )()( )( 11    向前差商 向后差商 ),(),()( 111111   nnnnnnn yxhfyyyxfxy 欧拉公式(隐式) • 2.积分法: )()()( 1 1 nn x x xyxydxxy n n    dxyxf n n x x ),( 1    ),(),(),( 11 nnnn x x x x yxhfdxyxfdxyxf n n n n    ),(1 nnnn yxhfyy   欧拉公式(显式) ),(),(),( 1111 11     nnnn x x x x yxhfdxyxfdxyxf n n n n ),( 111   nnnn yxhfyy 欧拉公式(隐式) )),(),(( 2 ),( 11 1   nnnn x x yxfyxf h dxyxf n n )),(),(( 2 111   nnnnnn yxfyxf h yy 梯形公式(隐式) 取初值: ),(1 nnnn yxhfyy  )),(),(( 2 111   nnnnnn yxfyxf h yy (预报值) (校正值) 改进欧拉公式 例 设      1)0( 2 y xyy (1)分别写出求解此微分方程的欧拉格式和改进 欧拉格式;(2)取步长h=0.1,求y(0.2)的近似值,并 与精确解y(x)=e-x2作比较. 解: (1)欧拉格式: )2(1 nnnn yxhyy  )21(1 nnn hxyy  改进欧拉格式: 整理为: 预报: )21(1 nnn hxyy  校正: )( 111   nnnnnn yxyxhyy 1)2.01( 0001  yxyy(2)欧拉格式: 98.0)2.01( 112  xyy 98.0)2.0( 2  yy 9608.0)2.0( 04.0  ey (仅一位有效数字) 改进欧拉格式: 1)2.01( 001  xyy 99.0)(1.0 110001  yxyxyy 9702.0)2.01( 112  xyy 9606.0)(1.0 221112  yxyxyy 9606.0)2.0( 2  yy (有三位有效数字) 3.单步法 • 单步法的一般形式是: ),,,( 11 hyyxhyy nnnnn    φ不含 yn+1 是显式,含 yn+1 是隐式. 欧拉格式: ),(1 nnnn yxhfyy  ),(),,,( 1 nnnnn yxfhyyx  ),(),,( yxfhyx  是显式 改进欧拉格式: ),(1 nnnn yxhfyy  )),(),(( 2 111   nnnnnn yxfyxf h yy )),(),(( 2 1 ),,,( 111   nnnnnnn yxfyxfhyyx ))),(,(),(( 2 1 nnnnnn yxhfyhxfyxf  ))),(,(),(( 2 1 ),,( yxhfyhxfyxfhyx  是显式 4.截断误差与精度 • 若假定 y(xn)=yn (第n步精确值等于近似值). 称: 为局部截断误差. • 若 则称相应方法有p阶精度. 111 )(   nnn yxyT )( 1 1    p n hOT • 若不作假定 y(xn)=yn . 称: 为整体截断误差. 111 )(   nnn yxye • 对固定点xn=x0+nh若 (或 )时有: 即 则称相应的 数值方法是收敛的. 0h n 0lim 1   n n e )(lim 11 0    nn h xyy )())(,( 11 1    pp nnn hohxyxT • 若 则称: 为局部截断误差的主项. 1 ))(,( p nn hxyx 例如欧拉格式: ),(1 nnnn yxhfyy  ),()( 11 nnnnn yxhfyxyT   ))(,()()( 1 nnnn xyxhfxyxy   )()()( 1 nnn xyhxyxy   由泰勒公式: 2 1 2 )( )()()()( h y hxyxyhxyxy nnnn   2 1 2 )( h y Tn    结论:欧拉格式有1阶精度,主项为 2 2 )( h xy n 例 证明:改进欧拉格式有2阶精度. )),(),(( 2 111   nnnnnn yxfyxf h yy ))),(,(),(( 2 nnnnnnn yxhfyhxfyxf h y  证明 由泰勒公式: )( 2 )( )()()( 32 1 hOh xy hxyxyxy nnnn    )( nn xyy  ))(,()( xyxfxy 微分方程 假定: )(),(),( ))(,( )( xyyxfyxf dx xyxdf xy yx  )))(,()(,( nnnn xyxhfxyhxf 由二元泰勒公式: )()())(,())(,())(,( 2 hOhxyxyxfhxyxfxyxf nnnynnxnn  ))),(,(),(( 2 1 nnnnnnnn yxhfyhxfyxf h yy  ))()(,( nnn xyhxyhxf  于是: )()( 2 )()( 3 2 hOxy h xyhxy nnn  )()()( 2 hOhxyxy nn  )())())(,())(,(()( 2 hOhxyxyxfxyxfxy nnnynnxn  )()( 3 111 hOyxyT nnn   即改进欧拉格式至少有2阶精度. )))(,()(,( nnnn xyxhfxyhxf 由二元泰勒公式: ))()(,( nnn xyhxyhxf  于是: )(),())()(,()(),(2),()( 2 xyyxfxyyxfxyyxfyxfxy yyyxyxx  )(),(),( ))(,( )( xyyxfyxf dx xyxdf xy yx  )(),( xyyxf y   hxyxyxfhxyxfxyxf nnnynnxnn )())(,())(,())(,(   )())(,(2))(,( 2 2 nnnxynnxx xyxyxfxyxf h  )())())((,( 32 hOxyxyxf nnnyy  )())())(,()(()()( 32 hOhxyxyxfxyhxyxy nnnynnn  ))),(,(),(( 2 1 nnnnnnnn yxhfyhxfyxf h yy   )( 2 )()( 2 nnn xy h xyhxy   111 )( nnn yxyT 即改进欧拉格式有2阶精度. )( 6 )( 2 )( )()()( 432 1 hOh xy h xy hxyxyxy nnnnn      )())())(,()(( 4 4 3 hOxyxyxfxy h nnnyn    )()())(,( 4 1 )( 12 1 43 hOhxyxyxfxy nnnyn  主项 9.3.R-K方法 1. r级与r阶R-K方法 • 据拉格郎日定理 只要对斜率 作出合理的估计. • 前述欧拉格式 • 前述改进欧拉格式 • 能否通过增多 来提高精度. ))(,()( )()( 1  yfy h xyxy kk   ))(,(  yf 1:  kk xx 其中 ),())(,( nn yxfyf  ),(),(( 2 1 ))(,( 11  nnnn yxfyxfyf  )(f    r i ii Kcyf 1 ))(,( 即 于是 2),,(),,( 1 1 1     iKhyhxfKyxfK i j jijnininn         r i i j jijnininn Khyhxfchyy 1 1 1 1 ),(  )(单步法       r i i j jijnininn Khyhxfchyx 1 1 1 ),(),,(  ).(,,, 显式称为待定常数ijiic  ).(, 显式称阶精度若达r .,, 格式达最高精度使选择 KRc ijii  r级R-K格式 r阶R-K格式 2.2阶R-K格式 )( 22111 KcKchyy nn 格式级 KR2 )),(,(),,( 21221 nnnnnn yxfhyhxfKyxfK   )()(),(),(),( 2 2122 hOxyyxfhyxfhyxfK nnnynnxnn   )())(),( ),(()()( 32 21 22211 hOhxyyxf yxfcxycchyy nnny nnxnnn     )( 2 )( )()()( 32 1 hOh xy hxyxyxy nnnn    )())(),(),(( 2 1 )()( 32 hOhxyyxfyxfhxyxy nnnynnxnn  :,2 必须阶精度为达 ,121  cc , 2 1 22 c . 2 1 212 c 格式中满足该要求为级因此 KR2, , 2 1 21  cc当 ,12  即时的格式,121  )( 2 211 KK h yy nn  ),(),,( 121 hKyhxfKyxfK nnnn  .是改进欧拉格式 2阶R-K格式 ,1,0 21  cc当 2 即时的格式, 2 1 21  21 hKyy nn  ) 2 , 2 (),,( 121 K h y h xfKyxfK nnnn  .称为中点公式 3.3阶与4阶R-K格式 )( 3322111 KcKcKchyy nn  ),(1 nn yxfK  ),( 12122 KhyhxfK nn   ))(,( 23213133 KKhyhxfK nn   ,1321  ccc , 2 1 3322   cc . 6 1 3223 c :当满足要求 .称为 ,212   ,32313   , 3 12 33 2 22   cc 3阶R-K格式 )22( 6 43211 KKKK h yy nn  ),(1 nn yxfK  ) 2 , 2 ( 23 K h y h xfK nn  ),( 34 hKyhxfK nn  ) 2 , 2 ( 12 K h y h xfK nn  经典4阶R-K格式: 例 2.00 0)0( 1       x y yy 分别用欧拉方法(h=0.025)、改进欧拉方法 (h=0.05)和经典4阶R-K格式(h=0.1)解微分方程 并与精确解比较. 解: 欧拉格式: )1(1  nnn yhyy 改进欧拉格式: )1(1  nnn yhyy )2( 2 11   nnnn yy h yy 经典4阶R-K格式: )22( 6 43211 KKKK h yy nn  11  nyK 1) 2 ( 12  K h yK n 1) 2 ( 23  K h yK n 1)( 34  hKyK n 数值解: 025.0975.01  nn yy欧拉格式 xk .025 .05 .075 .1 .125 .15 .175 0.2 yk .025 .049 .073 .096 .119 .141 .162 .183 改进欧拉格式: 05.0025.0975.0 11   nnn yyy 05.095.01  nn yy xk .05 .1 .15 .2 yk .048750 .095123 .13924 .18120 经典4阶R-K格式: )1.02.02.01.0( 6 1 43211 KKKKyy nn  11  nyK 1)05.0( 12  KyK n 1)05.0( 23  KyK n 1)1.0( 34  KyK n xk 0.1 0.2 yk .09516250 .18126910 精确解: 1  x ey 比较: x 欧拉 改欧 R-K 精确y(x) 0.1 .096312 .095123 .09516250 .09516258 0.2 .183348 .181198 .18126910 .18126925 ah=0.1 y1=0. do i=1,2 xk1=1-y1 xk2=1-(y1+ah*.5*xk1) xk3=1-(y1+ah*.5*xk2) xk4=1-(y1+ah*xk3) y2=y1+(xk1+2.*xk2+2.*xk3+xk4)*ah/6 write(0,*) y2 y1=y2 enddo pause end 4. 变步长R-K格式 5 1 2 11 2 2)(         h Cyxy h nn 5 111)( hCyxy h nn   ∵跨了两步,是 两步误差之和. 于是 16 1 )( )( 2 11 11      h nn h nn yxy yxy )( 15 1 )(: 1 2 1 2 11 h n h n h nn yyyxy  向后误差估计 h n h n yy 1 2 1  令给定误差界ε, • 若△>ε,将步长h折半,直至△<ε. • 若△<ε,将步长h加倍,直至△>ε. 9.4单步法的收敛性与稳定性 1.局部截断误差与整体截断误差的关系 • 定理 设显式单步法 (1) 有p阶精度; (2) 关于y满足Lipschits条件 即 ; (3) 则:整体截断误差 . ),,(1 hyxhyy nnnn  ),,( hyx 2121 ),,(),,( yyLhyxhyx  00 )( yxy  )( p n hOe  :证 :证 )),(,()(1 hxyxhxyy nnnn 设 1 11 1 111 ),()(      p n p nnn hcThOyxyT 且 1111111 )()(   nnnnnnn yyyxyyxye ),,()),(,()(11 hyxhyhxyxhxyyy nnnnnnnn    )),,()),(,(()( hyxhxyxhyxy nnnnnn   ),,()),(,()(11 hyxhxyxhyxyyy nnnnnnnn    nnnn ehLyxyhLe )1()(  n p nnnn ehLhcyyTe )1( 1 11111    ))1()(1( 1 1 1 1 1    n pp ehLhchLhc 1 21 1 )1()2(    n p ehLhLhc ))1(()1()2( 2 1 1 21 1    n pp ehLhchLhLhc 2 321 1 )1())1()1(1(    n p ehLhLhLhc 1 11 1 )1())1()1(1(    kn kkp ehLhLhLhc  0 211 1 )1())1()1(1( ehLhLhLhc nnp    )0( 1)1( 0 2 1 1      e hL hL hc n p  )()1)1(( 21 ppn hOhhL L c   LhLnhL n   )2(~1)1( 2 .证毕 .0lim 1 收敛  n n e • 例如欧拉格式 ,若 关于y 满足Lipschits条件 关于y满足Lipschits 条件 收敛. • 而 关于y是很容易满足 条件的, 只要 连续即可. ),(),,( yxfhyx  ),( yxf ),,( hyx  ),( yxf Lipschits ),( yxf y 2.相容性 显式单步法 ),,(1 hyxhyy nnnn  假定 y(xn)=yn 111 )(   nnn yxyT ))),(,()(()( hxyxhxyhxy nnnn  )),(,()()( hxyxhxyhxy nnnn  ))()()0),(,(( 2 )( )( 2 0 2 hOh h xyxhh y hxy nnn          )())0),(,()(( 2 hOhxyxxy nnn   于是: 即 • 一般:若 就称该单步法 与初值问题相容. • 定理:显式单步法是 阶的. 它是相容的. 0)0),(,()()( 2 1  nnnn xyxxyhOT  )0),(,())(,()( 2 1 nnnnn xyxxyxfhOT  )0),(,())(,( xyxxyxf  1p  3.稳定性 • 以上假定 是离散格式的精确解.由于电脑的 字长有限,设求得的解由4舍5入得到,设为 , 其中 是误差.由于这个误差,会使以 后计算出的 产生误差 ,若 , 则称该方法是稳定的. • 稳定性很复杂,对某格式,常用模型方程 来研究. ny ny ~ nnnn yy  ~ )( nmym  m nm   yy  0 1)0( 100 : 100        xey y yy 精确解设例 ),(1 nnnn yxhfyy 欧拉法 nnn yhyhy )1001()100(  解: 若: ;0:01.0 稳定 nyh nn yyh  1:02.0 ;,01,1,01,1 3210 不稳定 yyyy nn yyh 2:03.0 1     ;,21,02,1 11110 不稳定   nn nyyy yy 一般: • 稳定性不仅与方法有关,也与h有关. • 一般: 相容+稳定→收敛 yy  nn yhy )1(: 1 欧拉法 nn yhy ~)1(~ 1  nnn yy ~ nn h  )1(1  .,11 格式稳定时当  h 9.5线性多步法 knknknnn yyyyy  121 ,,,,: 多步法 线性多步法        k i ini k i inikn fhyy 0 1 0  :)( 步法又称线性k .,, 是待定系数不全为零其中 ii  .,0,,0 为隐式多步法即显式多步法  kk  ),( ininin yxff   假定 yn+i= y(xn+i) i=0,1,2, …,k-1 • 误差 为局部截断误差. • 若 称k步法有p阶精度. • 若 称该数值解法与原微分方程是相容的. • 在 处作泰勒展开,可以选择系数 , 使计算格式有尽可能高的精度. knknn yxyT   )(1 )( 1 1    p n hOT 1p 0x ii  ,  )( 2 )( )()()( 2 nnnn xy ih xyihxyihxy )())(,(),( ihxyihxyihxfyxf nnninin   )( 2 )( )()( 2 nnn xy ih xyihxy 常用的线性多步法     k i iniknkn fhyyAdams 0 1:)(.1 格式阿当姆斯     3 0 4:)ln(.2 i ininn fhyyeMi 格式米尔尼 )4( 3 :)(.3 212   nnnnn fff h yySip 格式辛普生 )( :)min(.4 33221122 1103     nnnn nnn fffhya yayaygHam  格式汉明 5.预测-校正方法: 用显式给一个初始的,粗糙的近似值,称预测; 再用隐式一步迭代得到一个较精确的近似值, 称校正.注意:取同阶的显式格式与隐式格式相 匹配.