第八章基于MATLAB的科学计算—常微分方程数值解法

科学计算—理论、方法及其基于MATLAB的程序实现与分析微分方程（组）数值解法§1常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律：000022xtxxtxtFdtxdm（１）和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：00yytfaydtdy（２）的解为：dfeyetyttaat00（３）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：000,ytytttytFdtdyf（７）其中tytytytyn21（８）ytfytfytfytFn,,,,21（９）常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。§１.1常微分方程（组）的Cauch问题数值解法概论假设要求在点（时刻）kt，nk,,2,1处初值问题（７）的解00,,yttyy的（近似）值,如果已求得kt时刻的值kty或它的近似值ky(如0t时刻的值0ty)，那么将式（７）的两端在区间1,kktt上积分dttytftytydtdtdykkkkttkktt11,1(10)可得dttytftytykkttkk1,1(11)或dttytfyytykkttkkk1,11(12)显然,为了利用式（11）或（12）求得1kty的精确值（近似值1ky），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解,可以说,不同的式值积分方法将给出不同的Cauch问题的数值解法。§1.2最简单的数值解法——Euler方法假设要求在点（时刻）khttk0，ntthf0，nk,,2,1处初值问题（７）的解tyy的近似值。首先对式（７）的两端积分，得dttytftytydtdtdykkkkttkktt11,1（13）对于式（13）的右边，如果用积分下限ktt处的函数值kktytf,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面积），则有kkttkktythfdttytftytykk,,11（14）进而得到下式给出的递推算法—Euler方法001,,2,1,ytynkythfyykkkk（15）例１用Euler方法解如下初值问题，取3.0h，1030sin23yttyydtdy解：由（15）得1010,,2,1sin6.03.131yktyyykkkk结果如下：openEuler_Method.m00.511.522.530.70.80.911.11.21.31.4EulerMethodofSolvingInitialValueProblemIntegralCurveEulerCurve如果取1.0h，其结果如下图所示：Euler_Method00.511.522.530.70.80.911.11.21.31.4EulerMethodofSolvingInitialValueProblemIntegralCurveEulerCurve§1.3改进的Euler方法对于（15）的右边，如果被积函数用积分限ktt和1ktt处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有111,,2,1kkkkttkktytftytfhdttytftytykk（16）进而得到下式给出的递推算法：00111,,2,1,,2ytynkytfytfhyykkkkkk（17）通常算法（17）比Euler方法（15）的精度高，但是，按算法（17）求1ky时要解（非线性）方程（组），这是算法（17）不如Euler方法的方面，为了１）尽可能地保持算法（17）精度高的优点；２）尽可能地利用Euler方法计算简单的长处；人们采取了如下的称之为改进的Euler方法的折衷方案：预测kkkkythfyy,01（18）修正nkytfytfhyykkkkkk,,2,1,,20111（19）00yty例２Euler方法与改进的Euler方法的比较下图是当3.0h时比较的结果：openImproved_Euler_Method.m00.511.522.530.70.80.911.11.21.31.4CompareEulerMethodWithImprovedEulerMethodIntegralCurveEulerCurveimprovedEulerCurve§1.4Euler方法和改进的Euler方法的误差分析由Taylor公式22111,!21hOhtytftytttytttytytykkkkkkkkkkk（19）说明Euler方法的截断误差是2hO，类似地，由321!21,hOhtyhtytftytykkkkk（20）321111!21,hOhtyhtytftytykkkkk（21）以及hOtytyhtytytykkkkk11(22)让式（20）的两端减式（21）的两端，可得31132111,,241,,2hOtytftytfhtyhOhOhtytftytfhtytykkkkkkkkkkk（23）从上述推导Euler方法、改进的Euler方法的过程以及例１、例２容易看出，改进的Euler方法Euler方法的精度高，其原因在于：1在推导Euler方法时，我们是用待求解函数tyy在一点处的变化率kktytf,代替tyy在区间],[1kktt上的平均变化率：tythfdttytfkktt,,1（24）2而在推导改进的Euler方法时，我们是用待求解函数tyy在两点处变化率的平均值11,,21kkkktytftytf代替tyy在区间],[1kktt上的平均变化率；显然，通常11,,21kkkktytftytf比kktytf,更接近于tyy在区间],[1kktt上的平均变化率。由此启发人们：适当地选取区间],[1kktt上函数tyy若干点处的变化率，用它们加权平均值代替tyy在区间],[1kktt上的平均变化率，近似解的精度应更高。下面将要介绍的Runge—Kutta法就是基于上述想法得到的。§2Runge—Kutta法Runge—Kutta法是按选取区间],[1kktt上函数tyy变化率ytf,的个数的多少和截断误差的阶数mhO来区分的一系列方法，如１二阶的Runge—Kutta法（改进的Euler方法）21111212,,KKhyyhKytfKytfKkkkkkk（25）２三阶的Runge—Kutta法32112312146,2,2,KKKhyyhKyhtfKKhyhtfKytfKkkkkkkkk（26）３四阶的Runge—Kutta法1)古典形式432113423121226,2,22,2,KKKKhyyhKyhtfKKhyhtfKKhyhtfKytfKkkkkkkkkkk（27）2)Gill公式(具有减小舍入误差的优点)432113242131212222622222,222212,22,2,KKKKhyyhKhKyhtfKhKhKyhtfKKhyhtfKytfKkkkkkkkkkk（28）4Runge—Kutta法的一般形式niiikkkkkKahyhythyy11,,(29)其中hytkk,,称为增量函数(IncrementFunction)以及11,111,1122212123111121,,,,nnnnknknkkkkkkhKqhKqyhptfKhKqhKqyhptfKhKqyhptfKytfK（30）需要特别指出的是：在确定（29）、（30）中的参数ia，ip和ijq时，应该使（29）的右端和适当阶的（如n阶）Taylor展式一致，这样，至少对于底阶的R-K法来说，参与加权的斜率个数n与方法的阶数是一致的。例如：一阶的R-K法即Euler法，局部的截断误差（TruncationError）是2hO，所以，当微分方程的阶是一次函数时是精确的；二阶R-K法即改进的Euler法，局部的截断误差是3hO，所以，当微分方程的阶是二次函数时是精确的；类似地，四阶R-K法局部的截断误差是5hO。但是，一般五阶以及五阶以上的R-K法的局部截断误差不再具有上述的特点，因此用“性价比”来衡量，四阶R-K得到了广泛的应用。下面是Butcher'sFifth-OrderR-K方法（1964）：6543117321232790KKKKKhyykk（28）543216415324213121787127127273,169163,4321,28181,44,4,hKhKhKhKhKyhtfKhKhKyhtfKhKhKyhtfKhKhKyhtfKKhyhtfKytfKkkkkkkkkkkkk（29）§4线性多步法与常微分方程数值解法的分类一般的常微分方程初值问题的数值解法都是以递推的形式给出的，即是递推算法，根据递推算法，在计算1ky时，已经得到了前面各时刻的近似值iy，ki,,1,0，前面介绍的各种数值解法都有一个共同点：在计算1ky时，只用到了前一个时刻（当前时刻kt）的“信息”ky预测未来某时刻1kt系统的“状态”1kty，这样的数值解法称为单步法，对于一个动态过程tyy在1kt时刻的状态1kty而言，不仅前一个时刻的信息ky对它有影响，前若干个时刻的信息通常对它也有影响，显然，单步法的缺