1、第八章第八章 常微分方程数值解常微分方程数值解/* Numerical Methods for Ordinary Differential Equations */y 为人口数量，为人口数量，r为人口增长率，为人口增长率， y0为时刻为时刻x0的人口数量。的人口数量。 000)()(yxyxxryy 求常微分方程求常微分方程解析解解析解的方法有多种多样，但是利用这些方法，的方法有多种多样，但是利用这些方法，我们只能对极少数特殊类型的常微分方程求解；科学研究和我们只能对极少数特殊类型的常微分方程求解；科学研究和工程技术上的大量常微分方程的求解需借助于工程技术上的大量常微分方程的求解需借助于数值计
2、算方法数值计算方法。 比如，描述人口增长的著名比如，描述人口增长的著名人口模型：人口模型：常微分方程常微分方程是常用的数学模型。是常用的数学模型。该方程有该方程有解析解解析解y(x)=y0er(x-xo)。 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条条件件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x
3、) 和和 y2(x) 都成立，则上述都成立，则上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 要计算出解函数要计算出解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距hi=xi+1-xi称为称为步长步长，当，当 hi = h 为常数时称为为常数时称为等步长等步长。利用利用y0求节点求节点x1处的近似值处的近似值y1, 再从再从y1来求出来求出y2, 直至直至求出所有的求出所有的yn. 称之为称之为步进法步进法. 1 欧拉方法欧拉方法 /* Eulers Method *
4、/ 欧拉公式：欧拉公式：xixi+1向前差商近似导数向前差商近似导数hxyxyxyiii)()()(1 )1,., 0(),(1 niyxfhyyiiii对于对于xia,b, 有有 )(,()(iiixyxfxy )(,()()(1iiiixyxhfxyxy 则则 有有 将上式中的函数值将上式中的函数值y(xi)都用近似值都用近似值yi来表示来表示, 则有数值计算则有数值计算格式格式7.1 Eulers Method例例7.1 求解初值问题求解初值问题 取步长分别为取步长分别为h=0.1和和h=0.05, 进行计算进行计算, 结果有结果有 1010 ,2)(yxyxyy解解: 该方程的解析解是
5、该方程的解析解是y=(1+2x)1/2. 欧拉格式是欧拉格式是)2(1iiiiiyxyhyy k xi y(xi) yi (h=0.1) yi (h=0.05) 012345.0.00.10.20.30.40.5.1.01.10001.19181.27741.35821.4351.1.01.09771.18761.27131.35021.4174.1.01.09541.18321.26491.34161.4142.定义定义在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，考虑步计算是精确的前提下，考虑的截断误差的截断误差 Ri = y(xi+1) yi+1 称为称为局部
6、截断误差局部截断误差 /* local truncation error */。定义定义若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶阶精度。精度。 欧拉法的局部截断误差：欧拉法的局部截断误差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR )()(322hOxyih 欧拉法具有欧拉法具有 1 阶精度。阶精度。Ri 的的主项主项/* leading term */ 局部截断误差：局部截断误差：7.1 Eulers Method 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* impli
7、cit Euler method */向后差商近似导数向后差商近似导数hxyxyxyiii)()()(11 xixi+1y(xi+1)yi+hf(xi+1, y(xi+1)1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxyih 即隐式欧拉公式具
8、有即隐式欧拉公式具有 1 阶精度。阶精度。7.1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，不便于实际计算。公式，不便于实际计算。)()(311hOyxyRiii)1,., 0(),(1 niyxfhyyiiii)1,., 0(),(111 niy
9、xfhyyiiii)()(322hOxyRihi )()(322hOxyRihi 7.1 Eulers Method)(3hORi 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii7.1 Eulers Method 2/ )()
10、,(),(11cpipiiciiipyyyyxhfyyyxhfyy或者或者 例例7.2 求解初值问题求解初值问题 取步长分别为取步长分别为h=0.1和和h=0.05, 进行计算进行计算, 结果有结果有 1010 ,2)(yxyxyy解解: 该方程的解析解是该方程的解析解是y=(1+2x)1/2. 改进欧拉格式改进欧拉格式是是)2( ),22(211111iiiiiiiiiiiiiyxyhyyyxyyxyhyy k xi y(xi) yi(h=0.1) yi(h=0.05) 012345.0.00.10.20.30.40.5.1.01.09591.18411.26621.34321.4164.1
11、.01.09561.18351.26531.34211.4148.1.01.09541.18321.26491.34161.4142.7.1 Eulers Method2 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */建立高精度的单步递推格式。建立高精度的单步递推格式。hxyxyii)()(1)(hxyi )(,(hxyhxfii *)()(,()()(hKxyhxyhxhfxyxyiiiii 1xnxn+1xn+hK*)(,(iiiixyxhfyy1考察欧拉公式考察欧拉公式 K* K1=f(xi, yi)改进欧拉公式改进欧拉公式 ),(),(11212112hK
12、yxfKyxfKKKhyyiiiiiiK* (K1+K2)/2斜率斜率一定取一定取K1 K2 的的平均值平均值吗吗？步长一定是步长一定是一个一个h 吗吗？这里称为函这里称为函数数y(x)在区间在区间xn, xn+1上的上的平均斜率平均斜率. 几何意义是几何意义是y(x)在区间在区间xn, xn+1上点上点xn+h处的处的斜率斜率. )(,(*hxyhxfKii 根据微分中值定根据微分中值定理理, 存在存在01 2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设
13、下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入公式，得到代入公式，得到 )()()()()()()
14、()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 2 Runge-Kutta MethodStep 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库
15、库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKh
16、KyhxfKhKyhxfKyxfKKKKhyy 常用的三阶常用的三阶龙格龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),()(213121321122246hKhKyhxfKKhyhxfKyxfKKKKhyyiiiiiiii2 Runge-Kutta Method 常用的四阶常用的四阶龙格龙格-库塔法库塔法 ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 例例7.3 用四阶龙格库塔公式求解例题用四阶龙格库塔公式求解例题7.1 取取
17、h=0.1计算计算, 结果有结果有解解: 对于例题对于例题7.1有以下四阶龙格库塔计算格式有以下四阶龙格库塔计算格式k xi y(xi) yi K1 K2 K3 K4 334223112143211)(22222222)22(6hKyhxhKyKKhyhxKhyKKhyhxKhyKyxyKKKKKhyynnnnnnnnnnnnnn , ,012345.0.00.10.20.30.40.5.1.00.912870.845120.790570.74536.0.954760.878180.817470.767810.72622.1.01.09541.18321.26491.34161.4142.0.
18、952290.876050.815620.766200.72481.0.912620.844940.790390.745200.70697.1.01.09541.18321.26491.34161.4142.2 Runge-Kutta Method3 微分方程组与高阶方程微分方程组与高阶方程 /* Systems of Differential Equations and Higher-Order Equations */ 一阶微分方程组一阶微分方程组IVP的一般形式为：的一般形式为： )(,.),(,()(.)(,.),(,()(1111xyxyxfxyxyxyxfxymmmm初值初值000
19、2020101)(,.,)(,)(mmyxyyxyyxy 将问题记作将问题记作向量形式向量形式，令：，令： 001011.,.,.mmmyyyfffyyy 00)(),()(yxyyxfxy前述所有公式皆适前述所有公式皆适用于向量形式。用于向量形式。 改进欧拉公式的形式改进欧拉公式的形式 3 Systems of DEs and Higher-Order Equations 四阶龙格库塔公式四阶龙格库塔公式 )z ,y,x(hgzz)z ,y,x(hfyynnnnnnnnnn11预报预报：2211111111)z ,y,x(g)z ,y,x(ghzz)z ,y,x(f)z ,y,x(fhyyn
20、nnnnnnnnnnnnnnn 校正校正： 2262264321143211LLLLhzzKKKKhyynnnn)hLz ,hKy,hx(gL)hLz ,hKy,hx(fK)Lhz ,Khy,hx(gL)Lhz ,Khy,hx(fK)Lhz ,Khy,hx(gL)Lhz ,Khy,hx(fK)z ,y,x(gL)z ,y,x(fKnnnnnnnnnnnnnnnnnnnnnnnn334334223223112112112222222222220000 z)x( z),z ,y,x(gzy)x(y),z ,y,x(fy例例 高阶微分方程高阶微分方程3 Systems of DEs and Higher-Order Equations 10)1(1000)1()()(,.,)(,)(),.,(nnnnaxyaxyaxyyyyxfy化作化作一阶微分方程组一阶微分方程组求解。求解。引入新变量引入新变量)1(21,., nnyyyyyy ),.,(.1121nnnnyyxfyyyyy初值条件为：初值条件为：10102001)(.)()( nnaxyaxyaxy
