Chapter5_4_方程组

§6常微分方程组数值解法/*Linearmultistepmethod*/二、刚性(stiff)常微分方程组一、常微分方程组数值解法三、高阶微分方程设一阶方程组:2001212210012111)(),,(dd)(),,(ddytyyytftyytyyytfty如果用经典的R-K四阶方法求解,算式为:一、常微分方程组数值解法)22(6114131211111kkkkyynn)22(6124232221212kkkkyynn),,(21111nnnyythfk),,(21221nnnyythfk)2,2,2(212111112kykyhthfknnn)2,2,2(212111222kykyhthfknnn)2,2,2(222121113kykyhthfknnn)2,2,2(222121223kykyhthfknnn),,(232131114kykyhthfknnn),,(232131221kykyhthfknnn,2,1,0n一般地将两个方程的情形推广到m个方程的情形),...,2,1()(),..,,,(dd......),..,,,(dd),..,,,(dd002121222111miytyyyytftyyyytftyyyytftyiimmmmm记T21),...,,(myyyYT21),...,,(mfffFT020100),...,,(myyyY则方程的向量形式为：a≤t≤b00)(),(ddYtYytFtY于是经典的四阶R-K方法的向量形式为：,...)2,1,0n（),()2,2()2,2(),()22(61342312143211KYhthFKKYhthFKKYhthFKYthFKKKKKYYnnnnnnnnnn其中：mmnnnyyyY...21mnnmnnnnnnkkkYthfYthfYthfYthFK11211211),(),(),(),(mnnmnnnnnnkkkKYhthfKYhthfKYhthfKYhthFK222211121112)2,2()2,2()2,2()2,2(mnnmnnnnnnkkkKYhthfKYhthfKYhthfKYhtFK332312222123)2,2()2,2()2,2()2,2(mnnmnnnnnnkkkKYhthfKYhthfKYhthfKYhtFK442413323134),(),(),(),(其分量形式为：),...,,,,()2,...,2,2,2()2,...,2,2,2(),...,,()2(6132321314222212131212111221143211mmnnnniimmnnnniimmnnnniimnnnniiiiiiininkykykyhthfkkykykyhthfkkykykyhthfkyyythfkkkkkyy),...2,1(mi,...)2,1,0(n二、刚性(stiff)常微分方程组设常系数线性微分方程组)(ddtAYtY其中mmRAT21),...,,(myyyYT21))(),...,(),(()(ttttΦm假设A有m个互不相同的特征值λj和相应的特征向量rj(j=1,2,…,m),则通解为:mjjtjjtΨectYr1)()(假设A的所有特征值都有非负实部,即:0)(Rej且|)(Re|...|)(Re||)(Re|m21因此mjjtjjtrec1)(0mjjtjjrec1称为瞬态解.rect11称为快瞬态解.rectmm称为慢瞬态解.|)(Re)(Re|1mS称为微分方程组的刚性比.当S1时,即认为是刚性问题.用显式R-K法求解刚性问题,只有当步长很小时,才能保证其稳定性.应用隐式的R-K方法常收到好的效果.一个化学反应系统中提出的刚性问题的例子A，B，C是三个化学样本，Robertson反应如下：)()()(410710304.0fastCACBfastveryCBBBslowBA我们通过建模可以得到如下方程组324111004.0yyyy1)0(1y227324121031004.0yyyyy0)0(2y2273103yy0)0(3y令则方程组可写为：TyyyY),,(321)(YfYTY)0,0,1()0(对此Robertson系统而言，Jacobian矩阵为01060101061004.0101004.0272427342434yyyyyy在平衡点（1，0.0000365，0），Jacobian矩阵为021900365.0219004.0365.0004.0特征值为01405.026.21893显然当t→+∞时解的各个分量yi(t)(i=1,2,3)是指数衰减的,并趋于稳态解(y1,y2,y3)=(0,0,0),y1(t),y2(t),y3(t)趋于稳态解(0,0,0)的速度是由因子e-0.1t决定的.上述初值问题的精确解是:ttttteetyetyeety30000503502501.01)()()(假如试图利用四级Runge-Kutta方法求解上述初值问题,要求计算直至得到符合精度要求的稳态解为止.我们讨论计算过程可能遇到的问题:一.稳定性要求:二.为使解充分接近稳态解只需要:78.2ih3,2,1i3000034103000078.2h所以01.0te4041.0teet4041.0teet而实际上1t后5030000,ttee已经不起作用了!!!5440404041010tNh往后的计算我们当然希望使用大步长!但由于稳定性要求,仍要用小步长.从而耗费了巨大的计算量,并且误差积累的影响也随着计算步数的增加越来越严重.N是计算步数410N上述例子中,系数矩阵的特征值虽然都是负数,但绝对值相差非常悬殊.考虑n维非线性常微分方程组设是定义在[0,T]上的解,并满足0)0(YY)(tY现用表示(1)在附近的解,)()()(tZtYtY)(tY则应满足)(tZnnnRRTFRTtYTttYtFtY],0[:,],0[:)(0)),(,()((1)TttZtYtFtZY0),())(,()((2)称为(1)在处的线性化方程)(tY记矩阵的特征值为))(,()(tYtFtJY),,2,1(),(niti称解为局部稳定,否则是不稳定的.)(tY若Ttniti0,,,2,1,0)(Re(3)定义1设是方程组(1)的一个解.假定相应Jacobi矩阵J(t)的特征值满足(3),并且],0[)(Re)(ReTttmintmaxiiii则称(1)在附近为刚性方程组.)(tY)(Re)(Retmintmaxiiii为刚性比.刚性方程组的解快变部分慢变部分一般地说,隐型方法比显型方法具有更大的绝对稳定区域,因此使用隐型方法求解刚性方程组更为合适.隐型Runge-Kutta方法Adams内插方法设高阶微分方程为令)()(1xyxy)()(2xyxy)()()1(xyxynn……),,,,()()1()(nnyyyxfxybxa)1(0)1()(nnyay…0)(yay0)(yay则)()(21xyxy)()(32xyxy))(,),(),(,()(21xyxyxyxfxynn01)(yay02)(yay…)1(0)(nnyay三、高阶微分方程§7常微分方程边值问题的数值解法/*Linearmultistepmethod*/设二阶线性常微分方程为常见边界条件有三类：],[),()()(baxxryxqyxpy101011)()(,)()()(,)()(,)(bybyayaybyaybyay差分方程的建立对区间[a,b]作等距分划：),...2,1,0(njjhaxjnabh由数值微分公式)(62)()()()3(211jjjjyhhxyxyxy)(12)()(2)()()4(2211jjjjjyhhxyxyxyxy),(,11jjjjxx代入微分方程得差分方程：jjjjjjjjjryqhyyphyyy2211211)1,...2,1(nj这是求yj(j=0,1,…,n)的n-1个方程,还缺的两个方程由边界条件给出。第一类边界条件y0=α,yn=β11221112211122222211nnnnnnnnnncdddadyyyybacbacbacbjjjjjjjjrhdphcqhbpha22;212;21其中整理得:追赶法求解yj(j=1,…,n-1)第二、三类边界条件由数值微分公式将边界条件离散化：1121210234243hyyyhyyynnn第二类边界条件1012100210234243nnnnyhyyyyhyyy第三类边界条件将两个方程分别与差分方程组联立求解yj(j=0,1,…,n)