常微分方程组数值解

在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近似值。其中x是质量，m是离开平衡位o的距离，t为时间，c为弹簧系数。例如：弹簧一质量系统的振动问题经一定的简化后可用一个二阶常微分方程oxmcdtxd22来描述。在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一起组成定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。mxxoc:)()('00'0000成一个初值问题唯一确定。这就可以写律都确定时，弹簧振动规和初始速度时的初始位置当弹簧在振动开始时刻txxtxdtdxxtxtttt我们现在讨论常微分方程的数值解法。先从最简单的一阶常微分方程的初值问题出发开始讨论。)1()(),(),,(0yaybaxyxfdxdy由常微分方程理论可知：只要上式中的函数f(x,y)在区域G={a≤x≤b,－∞＜y＜∞}内连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L，使'00'0022)()(0xtxxtxxmcdtxd1,(0,1,1)iiihxxin两相邻点间的距离称为步长。)~1.(0niihxxhhii，有（常值）时称为等步长当至于初值问题（1）的数值解法，常采用差分方法，即把一个连续的初值问题离散化为一个差分方程来求解。具体地，将（1）离散化后，求找其解y=y(x)在一系列离散点bxxxxxani210niyyyyy210上的近似值下面分析均假定满足上述条件。2121),(),(yyLyxfyxf12,()Gyyyyx对内任意两个都成立，则方程(1)的解必定存在且唯一。对于初值问题（1），先将其离散化，即把[a,b]区间n等分，.),2,1,0(nabhniihaxi其中得各离散节点因为初值问题中的初始条件已知，即可利用已知的来求出下一节点处的近似值；再从来求如此继续，直到求出为止。这种用按节点的排列顺序一步一步地向前推进的方式求解的差分算法称为“步进式”或“递推式”算法，它是初值问题数值解法的各种差分格式的共同特点。因此，只要能写出由前几步已知信息来计算的递推公式（即差分格式），即可完全表达这种算法。0)(yay0y)(1xy1y1ynyiyyy101iy...2y若将和的近似值分别记为和，则得：（3）)(ixyiy)1,2,1,0(),(1niyxhfyyiiii))(,()()(1iiiixyxhfxyxy或)1,2,1,0().()()('1nixhyxyxyiii将上式改写成：无穷小量充分小时，可忽略高阶有界且当),(''2)(2iiyhhy展式为：处在为方程的解，则设Tayloryxxyxyyiii),()(1)2()1,2,1,0()(),(''2)(')()(121nixxyhxhyxyxyiiiiiii)(1ixy1iy这就是Euler公式（格式）。利用它可由初值出发逐步算出。这类形式的方法也称为差分方法。nyyy21iy)(iixyy1)(1ixyyi)1,1,0(),(),(2)(1211nixxyhyxyRiiiiiii)4()(2hO定义：如果局部截断误差为，则这种数值算法的精度为p阶，故Euler格式的精度为一阶。从几何意义上来看，如图，0y当假定为准确值，即在的前提下来估计误差，这种截断误差称为局部截断误差。ix)(1phO由（2）、（3）知Euler公式在处的局部截断误差为：由方程（1）知，其积分曲线y=f(x)上任一点(x,y)的切线斜率都等于函数f(x,y)的值。从初始点（即点）出发，作积分曲线y=y(x)在点上的切线（其斜率为）与直线相交与点（即点），得到作为的近似值，则有dxdy0p),(00yx10pp),(0oyxf),(11yx1y)(1xyynpOxnniixxxxxxx112100p1pip1ip1np'1p'2p'ip'1ip'1np'npy=f(x)0p1xx1p2p),()(0001010yxhfydxdyxxyyoxx10))(,()()(01xxdxxyxfxyxy与相比较知，这时用切线近似代替了曲线段点近似代替了点，近似代替了近似代替了。递推继续从点出发，作一斜率为的直线与直线相交于点（即点），得到作为的近似值。…如此直到点。这样得出一条折线近似代替积分曲线，当步数越多时，由于误差的积累，折线可能会越10pp1'10,ppp'1p1y),(),(001yxhfxy10))(,(xxdxxyxf1p),(11yxf21pp2p),(22yx2y)(2xynpnippppp210'''2'10nipppppnippppp2101xx曲线。来越偏离真解'''2'10nippppp1)0()10(2'yxyxyy例：求解初值问题解：为便于进行比较，我们后面将用多种数值方法求解上述初值问题。这里先用Euler公式，此处具体格式为：取步长为h=0.1，计算结果略。由结果可见Euler方法的精度很差。1.0,1)0(,0)2(001hyyxyxyhyynnnnn即为Euler格式（3）。)5(),2,1,0())(,()('nixyxfxyiii)1,2,1,0()()()('1nixyhxyxyiii用向前微商)1,2,1,0()),(,()()(51nixyxhfxyxyiiii）可得：代入（则得分别用其近似值代替，)(),(1iixyxy)1,2,1,0(),(1niyxhfyyiiiiix因为差商是微分的近似，所以Euler格式也可用差商近似代替导数的离散方法来得到。在节点处有：显然Euler格式具有递推性，在计算时只要用到前一步所得的结果一个信息就够了，因此是一种单步格式或称一步格式。若用不同的数值微分计算方法也可导出其它形式的算法。例如：用向后差商表示的数值微分公式1iyiy)1,2,1,0()(2)]()([1)(11'niyhxyxyhxyiiii))(21)())(()()((2111'1iiiiiiiixxyxxxyxyxyTaylor展式由可得：代入))(,()(111'iiixyxfxy12111),(2))(,()()(iiiiiiiixxfhxyxhfxyxy和局部截断误差：则可得近似计算公式：)6()1,2,1,0(),(111niyxhfyyiiii)7()1,2,1,0()()(),(2212nihoxxyhRiiiii（6）称为向后Euler公式，又称为隐式Euler公式(后退Euler格式)。后退Euler公式与Euler公式有着本质的区别，后者是关于的一个直接计算公式，这类公式称作显式的；而前者，即（6）中右端含有未知的，它实际上是关于的一个函数方程。这类公式称作隐式的。1ny1ny1ny显式与隐式两类方法各有特点，使用显式算法远比隐式算法方便，但考虑数值稳定性等因素，人们常选用隐式算法。隐试算法（6）常用迭代法来实现，而迭代过程实质上是逐步显式化。设用显式Euler格式算出作迭代初值，以此代入（6）右端，使之转化为显示，直接计算得：，再用代入（6）右端又有：),()0(1iiiiyxhfyy)0(1iy(0)(0)111(,)iiiiyyhfxy)1(1iy),()1(11)2(1iiiiyxhfyy若迭代过程收敛，则极限值必为隐式方程（6）的解，从而可获得后退Euler方法的解。如此反复进行可得序列。}{)(1hiy,2,1,0),,()(11)1(1kyxhfyykiiiki1)(1limikikyy)8()],(),([2111iiiiiiyxfyxfhyy从几何上看，梯形公式是取区间两端点处斜率的平均斜率。],[1iixx)(2iyh比较Euler公式和后退欧拉公式的截断误差公式（4），（7）易见两者为同阶，但符号相反，因此就想到将这两种算法进行算术平均，其结果可能消除误差的主要部分，而获得更高精度的结果。这种平均方法称为梯形方法，公式为：1ixxEuler方法是过点以斜率引直线交的点A。后退Euler方法是以点处的斜率为斜率，从点引直线交于另一点B。A、B两点都偏离点点，梯形方法就是取A、B两点的中点P作为Q2近似，上图表明梯形方法确实改善了精度。),(iiyxf),(112iiyxQ),(11iiyxf1Q1ixx2Q1iixxy=f(x)xyoAB2Q1QP),(1iiyxQ)2,1,0())(,()('nixyxfxyiii代入)(3))(,(2)()('311iiiiiyhxyxhfxyxy得从而也可得二步Euler公式及其截断误差为：)10()1,1,0()()(3)9()1,1,0(),(23'311nihoyhRniyxhfyyiiiiii),2,1,0()(6)]()([21)('211'niyhxyxyhxyiiii也可以由导数的中心差分近似式得到：将区间[a,b]n等分，得子区间在第i+1个子区间上，积分1,1,0],,[1nixxii11),()()()('1iiiixxixxiidxyxfdxxyxyxyEuler公式、后退Euler公式以及梯形公式，二步Euler公式均可利用前面的数值积分技术得到。axyxybaxyxfdxdy000)(],[),,(例如：对于初值问题：对右端利用左矩形公式可得即Euler格式))(,()()(1iiiixyxhfxyxy),(1iiiiyxhfyy式可得：上进行，利用中矩形公间如果将积分改成在子区],[11iixx格式。二步Euleryxhfyyiiii),(211格式隐式即可得：对右端利用右矩形公式Euleryxhfyyxyxhfxyxyiiiiiiii),()).(,()()(111111梯形公式式可得：对右端利用梯形积分公)],(),([2111iiiiiiyxfyxfhyy各公式的截断误差可直接利用数值积分截断误差估计而得。从而可知梯形公式（8）的截断误差为：)11()()11,0()(1231'3honixxyhRiiiii梯形公式也是隐式的，可用迭代法求解，与后退Euler方法一样，仍用Euler方法提供迭代初值，其迭代格式为：为分析迭代过程的收敛性，将（12）与（8）相减得：（12）k=0，1，2，…..L为f(x,y)关于y的Lipschitz常数.如果选取h充分小使得则.时有这表明迭代过程（12）是收敛于（8）的解的。kiiiiikiiiiiyxfyxfhyyyxhfyy111101,,[2,12hlk11ikiyy011111111111122],,[2