计算方法第八章(常微分方程数值解).

第八章常微分方程初值问题数值解法本章主要研究常微分方程初值问题的数值求解：通常，假设函数f关于第二个变量满足李普希茨条件（L条件），即为存在常数L0，使得0()(,())()yxfxyxaxbyay1212(,)(,)fxyfxyLyy第一节一般概念1.1欧拉法及其简单改进0()(,())()yxfxyxaxbyay方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的解的近似值。012111NNnnnnnaxxxxxbhxxxx称为到的步长（通常取为常数h）1()()()|nnnxxyxyxyxh欧拉方法100(,),()nnnnyyhfxyyxy，记yn为y(xn）的近似计算值，有例子：(01)(0)1yyxy10(,)(1)1nnnnfxyyyyhyhyyxye精确解为下面我们分别取步长为0.1与0.01进行计算，计算结果显示在下面的图中。步长为0.1的计算结果。步长为0.01的计算结果0.010.990050.990.10.904840.904380.20.818730.817910.30.740820.73970.40.670320.668970.410.663650.662280.590.554330.552680.60.548810.547160.90.406570.404730.910.402520.400680.990.371580.3697310.367880.36603DOUBLEPRECISIONh,y(0：100)OPEN(20,FILE='OUTPUT.DAT',STATUS=UNKNOWN)h=1.0/100y(0)=1.0do10i=1,100y(i)=y(i-1)*(1.0-h)write(20,*)i*h,y(i)10continueEND我们可得精度更高的欧拉公式：11()()()2nnnyxyxyxh22()()()()()212nnnyxhyxhyyxhOhh11002(,),()nnnnyyhfxyyxy欧拉中点公式欧拉公式中我们利用了近似公式1()()()nnnyxyxyxh光这个近似产生的误差为1()()1()()()2nnnyxyxyxyhOhh利用利用中点公式求解微分方程时，有一个问题，就是计算时需要两个迭代初值！这样的算法称为二步法。前面的欧拉法称为单步法。对于这个问题，我们可以先用欧拉公式，通过给定的初值计算出的值，然后再利用这两个值（y0和y1）进行计算，直到计算出全部节点上的值。111(,,,)nnnnnyxxyy一般的单步法：一般的k步法：1(,,,,,,)nknnknnknkyxxyyy1.2欧拉方法的其他改进微分方程数值解的关键在于对导数的处理，可以用差分来近似导数，也可以通过积分，将导数项化掉。对于方程：首先，作出划分设已经求出第n个节点的函数值，在区间上对方程两边积分容易看出，要求第n+1个节点的函数值，关键在于选择适当的积分公式计算积分！()(,())yxfxyx0121NNaxxxxxb1[,]nnxxny11()()(,())nnxnnxyxyxfxyxdx（1）如选择下矩形公式，则得这正是前面的欧拉公式。1(,)nnnnyyfxyh111(,)nnnnyyfxyh111[(,)(,)]/2nnnnnnyyfxyfxyh（2）如选择上矩形公式，则得这是所谓的后退欧拉公式。（3）如选择梯形公式，则得这是所谓的欧拉梯形公式。直接利用已经求得的已知节点上的值计算未知节点上的函数值的算法称为显式法。例如：欧拉公式、欧拉中点公式计算未知节点上的函数值时，用到了未知节点上的函数值，这种算法称为隐式法。例如：后退欧拉法、欧拉梯形公式显然，利用隐式法求微分方程的数值解是，需要从表达式中反解未知节点上的函数值。1.3隐式法的具体计算：例如欧拉梯形公式用迭代法计算yn+1的值。（1）简单迭代收敛的条件：11111[(,)(,)]/2(,)/2(,)/2nnnnnnnnnnnyyfxyfxyhyfxyhfxyh(0)1(,)nnnnyyhfxy(1)()111(,)/2(,)/2kknnnnnnyyfxyhfxyh/21Lh（2）牛顿迭代若用简单迭代，而且只迭代一步，这样组成的一组计算公式称为预测－－校正公式。（迭代初值称为预测，迭代步称为校正）()()(1)()11111()11(,)/2(,)/21(,)/2kkkknnnnnnnnkynnyyfxyhfxyhyyfxyh(0)1ny(0)1(0)111(,)[(,)(,)]/2nnnnnnnnnnyyhfxyyyfxyfxyh预测－校正公式也称为改进的欧拉法，将上面的组合公式改写为：注意到，将上式进一步改写为：这是我们最终使用的计算格式。11[(,)(,(,))]/2nnnnnnnnyyfxyfxyhfxyh1121211()2(,)(,)nnnnnnyyKKKhfxyKhfxhyK1nnxxh例子：取步长为0.1计算，结果如图。(01)(0)1yyxy11212111()2(,)(,)()nnnnnnnnyyKKKhfxyhyKhfxhyKhyK图：DOUBLEPRECISIONh,y(0:10),ak1,ak2OPEN(20,FILE='OUTPUT1.DAT',STATUS=UNKNOWN)h=1.0/10y(0)=1.0do10i=1,10ak1=-h*y(i-1)ak2=-h*(y(i-1)+ak1)y(i)=y(i-1)+(ak1+ak2)/2.010continuedo20i=0,10write(20,*)i*h,y(i),exp(-i*h)20continueEND同理，对于后退欧拉公式有预测－校正公式或改写为：111(,)nnnnyyfxyh(0)1(0)111(,)(,)nnnnnnnnyyhfxyyyfxyh11(,)(,)nnnnnnKhfxyyyfxyKh用此法解前面的例子步长0.1步长0.011.4误差估计定义：利用第n个节点或之前更多节点的函数精确值，利用近似公式数值计算第n+1个节点的近似值，所引起的误差，称为第n+1个节点上的局部截断误差。我们记为第n+1个节点上解的精确值，为假设等条件下计算所得的近似值，则局部截断误差为：如局部截断误差为，称为具有p阶局部截断误差。1ny1()nyx111()nnnyxy()pOh()nnyyx欧拉方法的误差分析：12()()()()()()()/2()()()/2nnnnnnnnyxyxyxhyxhhyxyxhyhyxyxyhh21()(,())()()(,())()nnnnnnnyxfxyxyxyxhfxyxOh1()(,())nnnnyyxhfxyx2111()()nnnyxyOh而完全类似的可以得到后退欧拉公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉梯形公式的局部截断误差为：2111()()nnnyxyOh3111()()nnnyxyOh3111()()nnnyxyOh定义：由初值计算到第n个节点的近似值与其精确值之间的误差称为第n个节点整体误差。定理：设下面求解微分方程的数值计算方法局部截断误差为p+1阶，且函数关于y满足利普希茨条件，同时初值是准确的，则整体截断误差为p阶。欧拉公式、后退欧拉公式的整体误差为1阶。欧拉中点公式、欧拉梯形公式的整体误差为2阶。1(,,)nnnnyyhxyh(,,)xyh(,,)(,,)xyhxyhLyy微分方程数值解法的进一步改进。再回到恒等式如果取作为节点，将被积函数用插值多项式来近似，用插值多项式带到积分中去求出积分，则可以得到所谓的亚当斯(Adams)显式公式局部截断误差：11()()(,())iixiixyxyxfxyxdx1,,,iiikxxx1011()iiiikikhyybfbfbfA1(1)[]()kkkiRyBhy类似地，如果取作为节点，可得亚当斯(Adams)隐式公式局部截断误差：11,,,,iiiikxxxx****10112()iiiiikikhyybfbfbfbfA*2(2)1[]()kkkiRyBhy进一步，如果我们将恒等式中的积分区间改为，并在此区间上用辛普森公式，1[,]iiixx111111111()()(,())()[(,())4(,())(,())]3iixiixiiiiiiiyxyxfxyxdxhyxfxyxfxyxfxyx1111[4]3iiiiihyyfff可得辛普森公式1.5绝对稳定性一个常微分方程数值解法应用于方程时对给定步长h称为绝对稳定，是指在某一步（如第一步）产生的误差（如计算机的存储误差），在计算中会逐步减小。称方程为试验方程，设在计算开始时产生误差（存储误差），此误差在以后会逐步减小，我们称该算法相对于是绝对稳定的，这样的的全体称为该算法的绝对稳定域。通常在复平面上考虑绝对稳定域。稳定域与实轴的交集为稳定区间yyhh(Re()0)yyh欧拉法的绝对稳定区域后退欧拉法的绝对稳定区域11h11h绝对稳定域越大，h的选取范围就越大，算法就越好！1.6局部截断误差的实用估计（1）用两种阶数相同的算法求解，计算出n+1步的近似值，从而得到局部截断误差估计。（2）用同样的公式，用不同步长计算出n+1步的近似值，从而得到局部截断误差估计。1.7隐式法隐式法具有较好的绝对稳定性！只不过在使用隐式法的时候，需要进行迭代，或者使用预测－校正计算格式。第二节泰勒级数法与龙格－库塔法对于方程：取计算步长为h，则，将函数进行泰勒展开如函数y(x)有p+1阶导数，容易得到p阶泰勒级数展开法：0()(,())()yxfxyxaxbyay1nnxxh21()()()()()/2!nnnnnyxyxhyxhyxhyx()21(1)12!!()(,)(1)!ppnnnnnppnyyyyhyhhpyExhhp公式中的导数用下面公式计算：例子：2(,)(,)(,)():(,)2(,)(,)(,)nnnnxnnynnnnnxxnnxynnnyynnynnnyfxyyfxyfxyyyfxyfxyyfxyyfxyy()cos04(0)1yxxxycos,sin,cosnnnnnnyxyxyx231cossin/2cos/6nnnnnyyhxhxhx步长0.1步长0.01龙格－库塔法：对于常微分方程的数值解法，一个关键在于选择精度高的算法计算下面公式中的积分。要高精度的计算积分，常用的方法是适当增加计算节点，不妨设用m个节点计算上面积分，节点为则积分为11()()(,())nnxnnxyxyxfxyxdx*,1,2,,inixxhim11(,())(,())nnmxininixifxyxdxhfxhyxh将积分改写为：则得公式：取这样的公式称为显式龙格