Ch3常微分方程的差分方法

1、第第六六章章 常微分方程的差分方法常微分方程的差分方法 考虑考虑一阶一阶常微分方程的常微分方程的初值问题初值问题: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条条件件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述初值问都成立，则上述初值问题的连续可微解题的连续可微解y(x)在在a, b 上上存在且唯一存在且唯一。| ),(),(|2121yyLyxfyxf 要计算出解函数要计算出解

2、函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距 为为步长步长，通常采用，通常采用等距节点等距节点，即取即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii推导公式推导公式: :(1) Taylor(1) Taylor展开法展开法将将)(1 kxy在在kxx 点进行点进行TaylorTaylor展开展开,! 2)()(,()()(121 kkkkkkkkkkxxhyxyxfhxyxy 忽略忽略2kh这一高阶项这一高阶项,)(kxy 由于由于因此因此分别用分别用,1 kkyy),(

3、kkkyxff 近似近似),(),(1 kkxyxy),(,(kkxyxf得得,1kkkkfhyy 结合初值条件结合初值条件00()y xy 即得即得(2.1).)2 . 2(00 () , yy x 1(,),kkkkkyyh f xy 1, 1 , 0 nk)1.2( )(,()(kkkxyxfxy ),()(yxfdxdyxy (2) 向前差商近似微分法向前差商近似微分法向前差商向前差商kkkhxyxy)()(1 近似近似),(kxy 得得)(,()()(1kkkkkxyxfhxyxy )3.2(将近似号改为等号，将近似号改为等号，),(,1kkkkkyxffyy 用用近似近似),(),

4、(1 kkxyxy),(,(kkxyxf并结合初始条件即得并结合初始条件即得(2.1)。 6.1 欧拉方法欧拉方法(3) 左矩形数值积分法左矩形数值积分法将将(1.1)两端从两端从kx到到1 kx积分积分,得得)1 . 1(),(yxfdxdy 1)(,()()(1kkxxkkdxxyxfxyxy数值积分采用左矩形公式数值积分采用左矩形公式 , 即即)(,()(,(1kkkxxxyxfhdxxyxfkk 由初始条件亦得由初始条件亦得(2.1).oxy)(,(xyxf)(,(kkxyxfkx1 kxkh用用),(,1kkkkkyxffyy 近似近似),(kxy),(,(kkxyxf得得,1kkk

5、kfhyy ),(1 kxy00 (), yfx ,1kkkkfhyy 1, 1 , 0 nk)1.2( 显式欧拉格式的几何意义：显式欧拉格式的几何意义：x0 x1)1,., 0(),(1 niyxfhyyiiii 在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提下，步计算是精确的前提下，考虑的截断误差考虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断误差局部截断误差。定义定义3.1.1: 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法，则称该算法有有p 阶精度阶精度。定义定义3.1.2: 欧拉法的局部截断误差：欧拉法的局部截

6、断误差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR 2322()()()hiyxO hO h 欧拉法具有欧拉法具有 1 阶精度。阶精度。称为称为显式欧拉格式显式欧拉格式,亦称为亦称为欧拉折线法欧拉折线法误差主项误差主项 如何如何改进改进？1. 隐式欧拉格式隐式欧拉格式向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称

7、为隐式欧拉格式隐式欧拉格式。简单情况可以把隐式格式化成显式格式；简单情况可以把隐式格式化成显式格式；一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。2. 隐式隐式欧拉法的局部截断误差欧拉法的局部截断误差11)(iiiyxyR)()(322hOxyih 即隐式欧拉格式具有即隐式欧拉格式具有 1 阶精度。阶精度。3. 梯形格式梯形格式 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注： 有局部截断误差有局部截断误差 ， 即梯形格式具有即梯形格式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法

8、有了进步。 但注意到该格式是但注意到该格式是隐式隐式格式，计算时不得不用到迭格式，计算时不得不用到迭代法，其迭代收敛性与欧拉格式相似。代法，其迭代收敛性与欧拉格式相似。)()(311hOyxyRiii4. 两步欧拉格式两步欧拉格式中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点格式具有即中点格式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2个初值个初值 y0和和 y1

9、来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为两步法两步法，而前面的三种算法都是而前面的三种算法都是单步法单步法 。6.2 预报预报-校正法校正法-改进欧拉法改进欧拉法Step 1: 先用先用显式显式欧拉格式作欧拉格式作预报预报，算出，算出),(1iiiiyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形格式的右边作梯形格式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预报预报-校正法校正法。可以证明该算法具有。可以证明该算法具有 2 阶精阶精度，同时可以看到它是个度，同时可以看到它是个单步单步递推格

10、式，比隐式格式的递推格式，比隐式格式的迭代求解过程迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高稳定性高于显式于显式欧拉法。欧拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii 例例 分别用显式分别用显式Euler方法，梯形方法和预报校正方法，梯形方法和预报校正Euler方法解方法解初值问题初值问题解：解：取取 h =0.1，1(,)kkkkyyhf xy (1)kkkyhyx 1),( xyyxf(1)kkh yhxh 91,1010010kky 0.1,00.1 .khxkhk , 1 xydxdy01x 1)0( y Euler

11、方法为方法为:梯形方法为：梯形方法为：111 (,)(,)2kkkkkkhyyf x yf xy 11(1) (1)2kkkkkhyyxyx 10.1,0.1 ,(1)0.1(1).kkhxkhkxkhk 10110521191 kyykk预报校正预报校正Euler方法：方法：1 . 00095. 0905. 01 kyykkkxky)(kkxyy )(kkxyy )(kkxyy kyky表表 6.1 数值例子表明，梯形方法和预报校正数值例子表明，梯形方法和预报校正EulerEuler方法比显式方法比显式EulerEuler方法有方法有更好的精度更好的精度。6.3 龙格龙格 - 库塔库塔(Ru

12、nge-Kutta)法法建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：1121211()2(,)(,)iiiiiiyyhKKKf xyKf xh yhK斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？首先希望能确定系数首

13、先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdx

14、dyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多

15、个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的