第8章 常微分方程数值解法 本章主要内容： 1．欧拉法、改进欧拉法 2

第8章 常微分方程数值解法 本章主要内容： 1．欧拉法、改进欧拉法 2 第8章 常微分方程数值解法 本章主要内容: 1(欧拉法、改进欧拉法. 2(龙格-库塔法。 3(单步法的收敛性与稳定性。 重点、难点 一、微分方程的数值解法 在工程技术或自然科学中，我们会遇到的许多微分方程的问题，而我们只能对其中具有较简单形式的微分方程才能够求出它们的精确解。对于大量的微分方程问题我们需要考虑求它们的满足一定精度的近似解的，称为微分方程的数值解法。本章我们主要 dy,f(x,y),,dx讨论常微分方程初值问题的数值解法。 , ,y(x)y,00, 数值解法的基本思想是:在常微分方程初值问题解的存在区间[a,b]内，取n+1个节点a=x,x,…,x=b (其中差h= x –x称为步长，一般取h为常数，即等步长)，在这01Nnnn-1 些节点上把常微分方程的初值问题离散化为差分方程的相应问题，再求出这些点的上的差分方程值作为相应的微分方程的近似值(满足精度要求)。 二、欧拉法与改进欧拉法 欧拉法与改进欧拉法是用数值积分方法对微分方程进行离散化的一种方法。 xn，1,y,f(x,y)，，y(x),y(x)，f(t,y(t))dt, 将常微分方程变为 n，1n,xn 1(欧拉法(欧拉折线法) 欧拉法是求解常微分方程初值问题的一种最简单的数值解法。 欧拉法的基本思想:用左矩阵公式计算(,)式右端积分，则得欧拉法的 ba,yyhfxynNh为: (,)(0,1,...,1),，,,,n，1nnnN 欧拉法局部截断误差 2h2,,R,y(,)x,,,x或简记为O(h)。 n，1n，1nn，1n，12 1 我们在计算时应注意欧拉法是一阶方法，计算误差较大。 欧拉法的几何意义:过点A(x，y)，A(x，y)，…，A (x ，y )，斜率分000111nnn 别为f(x，y)，f(x，y)，…，f(x ，y )所连接的一条折线，所以欧拉法亦称为0011nn 欧拉折线法。 例1用欧拉法解初值问题 dy,,,2xy(0,x,1),dx, ,y(0),1, ,0 (0.2) 1处的近似解。(计算过程保留4位小数)。 在x 【思路】 用欧拉法求解常微分方程的初值问题时，首先熟练掌握欧拉公式的一般形式， 根据具体题目写出找出欧拉公式的迭代式，并根据初始条件和所给步长进行迭代求解。 解 ? f(x，y),,2xy ，h,0.2， 欧拉公式为: y,y，hf(x,y),y，0.2(,2xy),(1,0.4x)yn，1nnnnnnnn (n,0,1,2,3,4,5) 列计算如下: n x y y(x) y(x)-y nnnnn 0 0 1 1 0 1 0.2 1 0.9608 -0.0392 2 0.4 0.92 0.8521 -0.0679 3 0.6 0.7728 0.6977 -0.0751 4 0.8 0.5873 0.5273 -0.06 5 1 0.3994 0.3679 -0.0315 2 2(改进欧拉法 改进欧拉法比欧拉法的计算准确，是对欧拉法的改进。改进欧拉法的基本思想:用梯形公式计算(,)式右端积分，则得改进欧拉法的计算公式为: hba, yyfxyfxynNh,,(,)(,)(0,1,...,1),，，,,,n，1nnnn，1n，1N2 利用改进欧拉法计算常微分方程初值问题时，我们应注意此公式为隐式表达式，需要对它进行迭代求解。计算时可以采用一次迭代和多次迭代，因此，就有改进欧拉法预估-校正法公式和反复迭代的改进欧拉法预估-校正法公式。 改进欧拉法预估-校正法公式: ，，0,y,y，hf(x,y)nnnn，1, ,h，，0,,,y,y，f(x,y)，f(x,y)(n,0,1,...,N,1)n，1nnnn，1n，12, 反复迭代的改进欧拉法预估-校正法公式: ，，0,y,y，hf(x,y)nnn，1n, ,h，，，，mm，1,,,y,y，f(x,y)，f(x,y)(n,0,？,N,1,m,0,1,...,)nnnn，1n，1，1n2, 3h3,,,R,y(,)x,,,x 改进欧拉法的局部截断误差或简记为O(h)。n，1n，1nn，1n，112 从局部截断误差的形式看，改进欧拉法是二阶方法，因此，它比欧拉法更精确。 例2用预估-校正法求初值问题 2,,y,,y,xy(0,x,1), , ,y(0),1, 在x=0(0.2)1的解。 【思路】掌握预估-校正法的计算公式，根据已知条件迭代求解。 2解 步长h=0.2，将代入预估-校正公式，整理得 f(x,y),,y,xy (0)2,y,0.8y,0.2xyn，1nnn, ,2(0)(0)2,y,0.9y,0.1xy,0.1(y，x(y)n，1nnnn，1n，1n，1, 列表计算如下: 3 y y(x ) n x nnn 0 0 1 1 1 0.2 0.80720 0.80463 2 0.4 0.63690 0.63145 3 0.6 0.49048 0.48918 4 0.8 0.37780 0.37720 5 1 0.29103 0.29100 (m)(m,1),3y,y,10 例3用改进欧拉法求解例1的初值问题,要求。 nn 【思路】掌握改进欧拉法的计算公式，根据已知条件迭代求解，并检验迭代解是否满足精度要求，若满足则确定此解为常微分方程在某点的近似解。 f(x,y),,2xy 解 将代入改进欧拉法的计算公式得: (0),y,y，hf(x,y),y，0.2(,2xy),(1,0.4x)yn，1nnnnnnnn, ,h，，mm,1，，，，m,1,,,,,y,y，f(x,y)，f(x,y),y,0.2xy，xy)nnnn，1n，1nnnn，1n，1n，12, 列表计算如下: ，，0，，1，，2，，3n ) y(x)-y y(xnnn x yyyy,ynnnnnn 0 0 1 1 1 1 1 0 1 0.2 1 0.96 0.9616 0.9615 0.9608 -0.0007 2 0.4 0.8846 0.8523 0.8549 0.8549 0.8521 -0.0028 3 0.6 0.7181 0.7003 0.7025 0.7022 0.6977 -0.0045 4 0.8 0.5337 0.5325 0.5327 0.5327 0.5273 -0.0054 5 1 0.3622 0.3750 0.3725 0.3730 0.3679 -0.0051 4 三、龙格-库塔法 1(龙格-库塔法 龙格-库塔法具有精度高、收敛、稳定，不需要计算高阶导数等优点，是求解微分方程初值问题的一组著名的显示单步方法，广泛应用于求解常微分方程的初值问题。本章我们介绍了二、三、四阶龙格-库塔法。 龙格-库塔法的基本思想: ,,yf(x,y)yxyx,(),()n，1n 在计算初值问题的数值解时，考虑均差，则由微分中值,,y(x)yh00, ,定理可得， 由初值问题可得公式为:y(x),y(x),hy(x，,h)(0,,,1)n，1nn 上式中称为区间y(x),y(x)，hf(x，,h,y(x，,h))f(x，,h,y(x，,h))k，1kkkkk上的平均斜率。如果给平均斜率一种计算方法，就可得到计算y(x )的近似值y的公n+1n+1式。 如果仅取处的斜率值作为平均斜率的近似值，则得到的的公式为欧拉yxf(x,y)nnnn，1公式; 如果取处的斜率值，的平均值作为平均斜率的近似x,xf(x,y)f(x,y)nn，1nnn，1n，1 值，则得到的的公式改进欧拉公式。 yn，1 ? 二阶龙格-库塔法的公式和局部截断误差: y,y，h(k，k),,,n，1n1122 ,,k,f(x,y) ,1nn , ,kf(x,yhk),，2n，1n1, 在上式中选择不同的参数，会得到不同的二阶龙格-库塔法公式，所以二阶龙格-库塔 3法公式不唯一。二阶龙格-库塔法公式的局部截断误差为Ο(h)。 常见的二阶龙格-库塔法公式有以下两种 改进欧拉法迭代公式 5 h,,,，，yy(kk)n，1n12,,2y,y，hk，n1n2,,,kf(x,y), k,f(x,y)1nn,1nn,,,，kf(x,yhk)hh2n，1n1,,k,f(x，,y，k)2nn122,, ? 三阶龙格-库塔法的公式和局部截断误差: 常见的三阶龙格-库塔法公式为 h,,，，，yy(k4kk)，n1n123,6, ,k,f(x,y)1nn, , hh,k,f(x，,y，k)2nn1,22, ,k,f(x，h,y,hk，2hk)3nn12, 4 三阶龙格-库塔法公式的局部截断误差为Ο(h)。 ? 四阶龙格-库塔法公式 ,, ,,yy,,yy，，hkhk，，nn11nn22,,kk,,ff((xx,,yy)),,11nnnn ,,hhhh,,,,，，，，kkff((xx,,yykk))22nnnn112222,, 通常所说的龙格-库塔法是指四阶龙格-库塔法，也称为标准龙格-库塔法。由于它是一步 法，(即已知，就可以求出，无需知道，，…的值)且它的计算精度yyyynn，1n,1n,2高，所以应用较多，但在计算时，因为每一步都需要计算四次f(x，y)的值，计算量较 大，所以，一般用来计算前几项的近似值，即“表头”。 四阶龙格-库塔法公式为的公式和局部截断误差: 6 h,y,y，(k，2k，2k，k)，n1n1234,6,k,f(x,y)1nn, ,hh, k,f(x，,y，k),2nn122, hh,k,f(x，,y，k)3nn2,22,kf(xh,yhk),，，,4nn3, 5 四阶龙格-库塔法的局部截断误差为Ο(h)。 四、单步法的收敛性和稳定性 1(收敛性 如果在无舍入误差且步长h充分小的情况下，求得的近似值足够精确地逼近真解 yn ，即:当时，一致地有 h,0y(x)y,y(x)nnn y,,y(x),y ?欧拉法整体截断误差:其中为真解，为在无舍入误差y(x)nnnnn情况下，从y用欧拉法计算公式求得的近似解。 0 ? 欧拉法的收敛条件:如果f(x,y)关于y满足Lipschitz条件，且局部截断误差R有 n 2h界，即则欧拉法收敛。且欧拉法的整体截断误差估计 R,M(n,1,2,...,N),n22 式为: hML(b,a)2,,(e,1) n2L ,,其中L为Lipschitz常数，b-a为求解区间的长度，。 Mmaxy(x),2a,x,b 3.稳定性和绝对稳定性 ?稳定性:指初始(或某步)产生的误差在后面的迭代计算中不会再扩大。即存在常数 ~~~,,,yyCyyC及h,0,h?h时，对任意两个初始值满足不等式 。 y,y00nn0000 ? 欧拉法稳定性的条件:如果f(x,y)关于y满足Lipschitz条件，则欧拉法稳定。 ~ ? 绝对稳定性:若对固定步长及任意两个初始值满足不等式 hy,y000 ~~,,,yyyy 。 nn00 7 ? 我们在讨论稳定性时应注意，一般在实际计算中只能取固定步长，它不可能任意缩小。所以绝对稳定性则表示的是对固定步长h，在初始(或某步)所产生的误差，在以0 后计算中不会逐步增长。由于绝对稳定性的成立和f(x,y)有关，讨论较为复杂。所以一般 ,yy,,,,地，我们对简单的微分方程的绝对稳定性进行讨论。使得数值方法绝对, ,y(x)y,00, 稳定的步长h和常数μ的取值范围称为绝对稳定域。(各种数值方法的绝对稳定域见课本)。 ,y,,20y,例4 对于初值问题， ,y(0),1, 证明当时，欧拉法绝对稳定。 h,0.1 证明 由欧拉公式得 y,,hy(120)nn,1 ~~y,,hy(120)nn,1 ne,1,20he,？,1,10he所以， ,10nn 1,20h,1,e,e当时，有 h,0.1n0 所以欧拉法绝对稳定。 本章学习以记忆公式和掌握证明题的推导方法为主。 8