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

第七章常微分方程的数值解法（二）第四节线性多步法第五节单步法的收敛性与稳定性第六节一阶方程组和高阶方程第七节边值问题的数值解法上页下页返回§4线性多步法用若干节点处的y及y’值的线性组合来近似y(xi+1).)...(...110111101kikiiikikiiiffffhyyyy--+---+++++++++=bbbbaaa其通式可写为：),(jjjyxff=当b-10时，为隐式公式;b-1=0则为显式公式.基于数值积分的构造法将在上积分，得到),(yxfy=],[1+iixx+=-+1))(,()()(1iixxiidxxyxfxyxy只要近似地算出右边的积分，则可通过近似y(xi+1).而选用不同近似式Ik，可得到不同的计算公式.+1))(,(iixxkdxxyxfIkiiIyy+=+1上页下页返回亚当姆斯显式公式利用k+1个节点上的被积函数值构造k阶牛顿后插多项式,有kiiifff--,...,,1]1,0[,)(+thtxNik+++=+1010)()())(,(1dthhtxRdthhtxNdxxyxfikxxikiiNewton插值余项dthtxNhyyikii)(101++=+/*显式计算公式*/局部截断误差为：+=-=++1011)()(dthtxRhyxyRikiii例1k=1时有)()(11--+=+=+iiiiiifftfftfhtxN--+-+=-++=10111)3(2)]([iiiiiiiiffhydtfftfhyydththtdxyfdhRxxi)1(!21))(,(1022+=)(1253iyh=上页下页返回注：一般有，其中Bk与yi+1计算公式中fi,…,fi-k各项的系数均可查表得到.)()2(2ikkkiyhBR++=10123k21231223245521-1216-2459-1252437249-12583720251fifi-1fi-2fi-3…Bk…………………Misprintonp.106常用的是k=3的4阶亚当姆斯显式公式)9375955(243211---+-+-+=iiiiiiffffhyy上页下页返回亚当姆斯隐式公式利用k+1个节点上的被积函数值fi+1,fi,…,fi-k+1构造k阶牛顿前插多项式。与显式多项式完全类似地可得到一系列隐式公式，并有，其中与fi+1,fi,…,fi-k+1的系数亦可查表得到.)()2(2ikkkiyhBR++=~kB~10123k21-21125249211282419121-245-241121-241-72019-fi+1fifi-1fi-2…Bk…………………~常用的是k=3的4阶亚当姆斯隐式公式)5199(242111--+++-++=iiiiiiffffhyy小于Bk较同阶显式稳定上页下页返回亚当姆斯预测-校正系统Step1:用Runge-Kutta法计算前k个初值；Step2:用Adams显式计算预测值；Step3:用同阶Adams隐式计算校正值.注意：三步所用公式的精度必须相同。通常用经典Runge-Kutta法配合4阶Adams公式.4阶Adams隐式公式的截断误差为)(72019)()5(511iiiyhyxy-=-++4阶Adams显式公式的截断误差为)(720251)()5(511iiiyhyxy=-++当h充分小时，可近似认为ii，则：19251)()(1111---++++iiiiyxyyxy)(270251)(1111++++-+iiiiyyyxy)(27019)(1111++++--iiiiyyyxyPredictedvaluepi+1Modifiedvaluemi+1Correctedvalueci+1Modifiedfinalvalueyi+1外推技术/*extrapolation*/上页下页返回.).,(),(2211尽可能高的阶，使该公式具有、展开原理确定试用其中给出一个三步显式格式例babaTayloryxffffhyynnnnnnn=++=--+上页下页返回§5收敛性和稳定性一、收敛性定义若某算法对于任意固定的x=xi=x0+ih，当h0(同时i)时有yiy(xi)，则称该算法是收敛的.例3就初值问题考察欧拉显式格式的收敛性.==0)0(yyyy解：该问题的精确解为xeyxy0)(=尤拉公式为iiiiyhyhyy)1(1+=+=+0)1(yhyii+=对任意固定的x=xi=ih，有iixhhxihyhyy])1[()1(/10/0+=+=ehhh=+/10)1(lim)(0ixxyeyi=上页下页返回.1)0(03*3的梯形法收敛证明：解初值问题例==+yyy上页下页返回0.00.10.20.30.40.5精确解改进尤拉法尤拉隐式尤拉显式节点xixey30-=二、稳定性例4考察初值问题在区间[0,0.5]上的解.分别用欧拉显、隐式格式和改进的欧拉格式计算数值解.=-=1)0()(30)(yxyxy1.0000-2.00004.0000-8.00001.6000101-3.20001011.00002.500010-16.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上页下页返回定义若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的.一般分析时为简单起见，只考虑试验方程yy=常数，可以是复数当步长取为h时，将某算法应用于上式，并假设只在初值产生误差，则若此误差以后逐步衰减，就称该算法相对于绝对稳定，的全体构成绝对稳定区域.我们称算法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考察隐式欧拉法11+++=iiiyhyyiiyhy-=[phone]++-=iih可见绝对稳定区域为：1|1|-h210ReImg注：一般来说，隐式尤拉法的绝对稳定性比同阶的显式法的好.上页下页返回例7隐式龙格-库塔法=++++=+++=+),...,1()...,(]...[11111mjhKhKyhxfKKKhyymmjjijijmmiibba而显式1~4阶方法的绝对稳定区域为++=+=+)2,2(1111KhyhxfKhKyyiiii其中2阶方法的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3---123ReImg无条件稳定上页下页返回§6微分方程组与高阶方程一、一阶微分方程组IVP的一般形式为：==))(,...),(,()(.........))(,...),(,()(1111xyxyxfxyxyxyxfxymmmm初值0002020101)(,...,)(,)(mmyxyyxyyxy===将问题记作向量形式，令：===001011...,...,...mmmyyyfffyyy==00)(),()(yxyyxfxy前述所有公式皆适用于向量形式.上页下页返回二、高阶微分方程====---10)1(1000)1()()(,...,)(,)(),...,,,(nnnnaxyaxyaxyyyyxfy化作一阶微分方程组求解.引入新变量)1(21,...,,-===nnyyyyyy===-),...,,(...1121nnnnyyxfyyyyy初值条件为：10102001)(...)()(-===nnaxyaxyaxy上页下页返回).2.0()2.0(2.0)(2)0(,1)0(1)(8yyxxyyyyxyyyxyKuttaRunge=====++--、处的函数值及其导数值在的解初值问题格式求二阶常微分方程用四阶例上页下页返回§7边值问题的数值解法2阶常微分方程边值问题===ba)(,)(),(),,(byaybaxyyxfy打靶法先猜测一个初始斜率y(a)=s，通过解初值问题===sayaayyyxfy)()(),,(y(b)=(s)找出s*使得(s*)=b，即把问题转化为求方程(s)-b=0的根.yx0abyx()b斜率=s0()s0斜率=s1()s1每计算一个(s)都必须解一个ODE.上页下页返回有限差分法将求解区间[a,b]等分为N份，取节点xi=a+ih(i=0,…,N)，在每一个节点处将y和y离散化.)(12)()()()()()4(2yhhhhxyxyhxyhxyxy-----+=泰勒展开)()()(2)(22hOhhxyxyhxy+-+-+=)(2)()()(2hOhhxyhxyxy+--+===-=-=+--+-+baNiiiiiiiyyNihyyyxfhyyy,1,...,1)2,,(2011211