常微分方程的数值解法

第六章常微分方程的数值解法§6.0引言§6.1算法构造的主要途径§6.2Runge-KuttaMethod算法§6.3线性多步法§6.4线性多步法的一般形式§6.5算法的稳定性、收敛性§6.0引言1．主要考虑如下的一阶常微分方程初值问题的求解:00,dyfxydxyxy微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。2.例如微分方程:xy'-2y=4x；初始条件:y(1)=-3。于是可得一阶常微分方程的初始问题24(1)3yyxy。显然函数y(x)=x2-4x满足以上条件,因而是该初始问题的微分方程的解。3.但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。因此，只能依赖于数值方法去获得微分方程的数值解。4．微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b]，初始点x0=a，将[a,b]进行划分得一系列节点x0,x1,...,xn，其中a=x0x1…xn=b。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk的近似值y(xk)，即y≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解。如图所示:§6.1算法构造的主要途径abx0x1x2...xn-1xn1欧拉公式1.1构造的思想：微分方程初值问题:00,dyfxydxyxy利用差商代替一阶导数，即010()()xxyxyxdydxh，则1000()()(,)yxyxfxyh。于是，可求出y(x1)的近似值y1，1000(,)yyhfxy同样地，可利用x1处的微分方程可得：2111(,)yyhfxy一般地，利用在xn处的微分方程可得：1(,),0,1,2,...nnnnyyhfxyn此式称为欧拉公式。1.2几何意义：对于微分方程y'=2(x+1),其通解是y=(x+1)2+c,是一个曲线族，当给定初值条件y(0)=2，其特解为y=(x+1)2+1。如图所示：由y(0)=2，过该曲线上一点(0,2)作曲线的切线，其斜率00,2xdyfdx，∴切线为：20,2(0)yfx,因此可计算出y1,如此，可根据：1111,()yyfxyxx，…故欧拉法又称欧拉折线法。1.3算例：例：2211dyydxxy解：h=0.2,xi=1+ihy1=y0+hf(x0，y0)=-1+0.2×2(1)211.0y2=y1+hf(x1，y1)=-1+0.2×2(1)20.93331.2y3=y2+hf(x2，y2)=-0.9333+0.2×2(0.9333)20.81.4y4=y3+hf(x3，y3)=0.8+0.2×2(0.8)20.61.6y5=y4+hf(x4，y4)=0.6+0.2×2(0.6)20.33331.8y6=y5+hf(x5，y5)=0.3333+0.2×2(0.3333)202精确解为：y=x2-2xxky(xk)ykek1.2-0.96-10.041.4-0.84-0.93330.09331.6-0.64-0.80.161.8-0.36-0.60.242.00-0.33330.33332.20.4400.44可以看出误差随着计算在积累。1.4Euler法的特点和误差迭代格式：1,0,1,2,,1nnnnyyhfxynN特点：(1)单步方法；(2)显式格式；(3)局部截断误差为2h。局部截断误差：当nnyyx时，由nyx按照欧拉方法计算来的1ky的误差称为局部截断误差。即，1[,]nnnnyxyxhfxyx是局部截断误差。如：2112nnnyxyxhyxhy欧拉法得：1,nnnnyyhfxy因此，局部截断误差是2h。2改进Euler法2.1方法构造在微分方程初值问题(,)dyfxydx,对其从kx到1kx进行定积分得:11()()(,())kkxkkxyxyxfxyxdx将右端的定积分用梯形公式来进行近似计算。111()()(,())(,())2kkkkkkhyxyxfxyxfxyx用1ky和ky来分别代替1()kyx和()kyx得计算格式:111((,)(,))2kkkkkkhyyfxyfxy这就是改进欧拉方法。2.2显式格式和隐式格式在欧拉式中1(,)kkkkyyhfxy每一步计算已知ky,直接用格式可以计算出1ky,此类格式称为显式格式。而在改进欧拉方法中111((,)(,))2kkkkkkhyyfxyfxy在每一步计算中1ky是未知,待求的,未知量在(,)fxy中这是一个方程,如f是非线形或超越函数,此方程是无法直接解出来(要依靠迭代法才能解出)。这类格式称为隐式格式。2.3算例例:(0)2yyxy用改进欧拉方法求解0.2h。解：111()2kkkkkkhyyyxyx解得:11122()1122kkkkhhyyxxhh注意：由于(,)fxyyx,是线形函数可以从隐式格式中解出1ky。问题的精确解是()1xyxexkx()kyxkyke欧拉方法误差0.22.4214032.4222220.0008130.021400.42.8918252.8938270.002000.051830.63.4221193.4257890.003670.09411……………2.010.3890610.438780.048721.19733预测－校正方法由于改进Euler法是隐式格式，无法从格式中直接求出1ny必须要解方程。下面用预测－校正方法来求隐式格式中的1ny。预测值:1,nnnnyyhfxy校正值:111,,2nnnnnnhyyfxyfxy此式相当于对隐式格式求1ky时采用迭代的方法,用欧拉格式得到的1ky作为初始值迭代公式迭代一次而已,此公式代入后得:1((,)(,(,)))2kkkkkkkkhyyfxyfxhyhfxy如改写成平均的形式为:11(,)(,)1()2pkkkckkpkpcyyhfxyyyhfxyyyy§6.2龙格-库塔法Runge-KuttaMethod1龙格-库塔法的思想1.1考虑微分方程的初值问题:00'(,)()yfxyyxy根据微分中值定理有：1()()'()kkkyxyxyxhh，其中0θ1。于是'()(,())kkkyxhfxhyxh即1(()()),)(kkkkfxyxyxhhyxh我们称(,())kkfxhyxh为y(x)在区间[xk,xk+1]上的平均斜率，记作K，其中，θ是存在但是未知的。因此，如何对平均斜率K进行近似计算，相应地就得到一种近似公式，或称为微分方程的一种计算格式。1.2例如：用f(xk,yk)作为平均斜率K的近似值就得到欧拉格式：1(,)kkkkyyhfxy。用111((,()(,))2kkkkfxyxfxy作为平均斜率K的近似值就得到比欧拉格式高一阶精度的格式，即，改进欧拉格式的预测-校正方法:111((,)(,))2kkkkkkhyyfxyfxy1.3启发(Motivative)能否在二维平面中x∈[xk,xk+1],y∈[yk,yk+1]上多找一些f(x,y)点，在这些点上作函数值的平均数并以此作为平均斜率K的近似值。由于自由度的增加，使得1()pOh的p能够提高，从而达到提高精度的目标，这就是龙格-库塔法的基本思想。1.4R-K的一般形式在11(,)|[,],[,]kkkkxyxxxyyy二维平面区域上取N个点，得到公式:11NkkjjjyyhK，其中Kj是二元函数f(x,y)在这N个点上的值。12221133311322,11,11(,),(,),(,()),(,())kkkkkkjkjkjjjjKfxyKfxhyhKKfxhyhKKKfxhyhKK其中,,jjj是待定常数。2二阶龙格-库塔法2.1计算格式:12221111122,,kkkkkkKfxyKfxhyhKyyhKK应用Taylor公式求参数：22112,,,。将f(x,y)在(xk,yk)上展开12221122211(,)()(,)(,)()kkkkkkkxyKfxyyxKfxhyhKfxyhfhKfOh又2311()()'()''()()2kkkkyxyxhyxhyxOh因'''()(,())'kkkxyxyyfdyxfxyxffyfffdx以上两式代入1()kyx中，又由112()kkyyhkk得23212221()()2(()xyxyhyhffffOhyhffhfhffOh比较两边得1222212:11:21:2xyffff四个未知量,三个方程,有无穷多组解，但局部截断误差均为3()Oh,是二阶精度的。例如:(1)取21,则211,1212,就是改进欧拉公式的预测-校正方法(或二阶龙格-库塔法)。121112,,()2kkkkkkKfxyKfxhyhKhyyKK(2)取212，则2112，则120,1,得二阶中点法：12112,11,22kkkkkkKfxyKfxhyhKyyhK3三阶龙格－库塔方法3.1计算格式:122211333113221112233(,)(,)(,())()kkkkkkkkKfxyKfxahyhKKfxahyhKKyyhKKK参数有12323,,,,,213132,,共有八个参数。将f和1kyx均在(,)kkxy点Tayor展开,有1,kkkKfxyyx2222223222122212111()22xyxxxyyyKfahfhffahfahffhffOh2233311322322333113223113221()212xyxxxyyyKfhfKKhffhfKKhfKKfh利用1Kf,2K按展开式代入(高于2h的不计)得:222333132232213222222233331323132()1122xyxyyxxxyyyKfhfhffhffhffhfhffhffh241126kkkkkhyxyxhyxyxyxh()kyxf,kkxydyxfxyxfffdx222kxyxxxyyyxyyddyxyxffffffffffffdxdx以上展开式代入1112233kkyyhKKK中比较得:f1231(1)xf223312(2)yff2123132312(3)xxf2222331112