第六章常微分方程数值解 / Numerical Methods for Ordinary Differential Equations * 6考虑一阶常微分方程的初值问题/ nitial-Value Problen*: d 女=f(x,y)xen,b y(a)=yo 只要∫(x,y)在{a,b×R1上连续,且关于y满足 Lipschitz条 件,即存在与x,y无关的常数L使|f(x,yn)-f(x,y2)|sL|-y2 对任意定义在[a,b上的y(x)和y2(x)都成立,则上述ⅣP存 在唯一解。 要计算出解函数y(x)在一系列节点a=x0x1.<xn=b 处的近似值v≈y(x;)(=1,…,n) 节点间距h=x+1-x(i=0…,n-1)为步长,通常采用等距节点, 即取h1=h(常数)
第六章 常微分方程数值解 /* Numerical Methods for Ordinary Differential Equations */ 考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: = = 0 ( ) ( , ) [ , ] y a y f x y x a b dx dy 只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使 对任意定义在 [a, b] 上的 y1 (x) 和 y2 (x) 都成立，则上述IVP存 在唯一解。 | ( , ) ( , )| | | 1 2 1 2 f x y − f x y L y − y 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 y y(x ) (i 1, ... ,n) i i = 节点间距 为步长，通常采用等距节点， 即取 hi = h (常数)。 ( 0, ... , 1) hi = xi+1 − xi i = n −
§1欧拉方注 k 亦称为欧拉折线法 向前去 / Euler's polygonal arc method*/ y(x)≈y((xn)=yn+h/(x1,)远为n0 y+1=y2+hf(x1,y)(i=0,…,n-1) 定义在假设y=y(x),即第i步计算是精确的前提下,考 虑的截断误差R1=y(x1)-y称为局部截断误差 F local truncation error*。 定义若其然 R;的主项 为Q(+),则称该算法有 阶精度 / leading term */ G欧拉法的局左: R=px)=1ye+小+号y(x)+O川-+订(x川 yx)+O()欧拉法具有阶精度
§1 欧拉方法 /* Euler’s Method */ ➢ 欧拉公式： x0 x1 向前差商近似导数 h y x y x y x ( ) ( ) ( ) 1 0 0 − ( ) ( ) ( ) ( , ) 1 0 0 0 0 0 y x y x + hy x = y + h f x y 1 y 记为 ( , ) ( 0, ... , 1) yi+1 = yi + h f xi yi i = n − 定义 在假设 yi = y(xi )，即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1 ) − yi+1 称为局部截断误差 /* local truncation error */。 定义 若某算法的局部截断误差为O(h p+1 )，则称该算法有p 阶精度。 欧拉法的局部截断误差： ( ) [ ( ) ( ) ( ) ( )] [ ( , )] 3 1 1 2 2 i i i i h i i i i i R = y x − y = y x + hy x + y x +O h − y + hf x y + + ( ) ( ) 3 2 2 y xi O h h = + 欧拉法具有 1 阶精度。 Ri 的主项 /* leading term */ 亦称为欧拉折线法 /* Euler’s polygonal arc method*/
欧拉公式的改进: a隐式欧拉法/ implicit Euler method 向后差商近似导数-y(x)yx)yx y(ps yo thf(xu,y(u Vi+i=y +h Hey! Isnt the leading term of the local truncation errorof Euler's method y(x ) 由于末为 Seems that we can make a good 倒到,故 称为隐式 use ofit 式/ explicit* 欧拉公式。 一般先用显式计算一个初值,再迭代求解。 隐式欧拉法的局部截断误差 R=y(x1)-y+=-2y(x1)+O(h3 即隐式欧拉公式具有1阶精度
➢ 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 h y x y x y x ( ) ( ) ( ) 1 0 1 − x0 x1 ( ) ( , ( )) 1 0 x1 y x1 y x y + h f ( , ) ( 0, ... , 1) 1 1 1 = + = − yi+ yi h f xi+ yi+ i n 由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。 隐式欧拉法的局部截断误差： 1 1 ( ) i = i+ − i+ R y x y ( ) ( ) 3 2 2 y xi O h h = − + 即隐式欧拉公式具有 1 阶精度。 Hey! Isn’t the leading term of the local truncation error of Euler’s method ? Seems that we can make a good use of it … ( ) 2 2 i h y x
平 梯形公式/ trapezoid formula*—显、隐式两种算法的平 Ji+FD;+If(x y)+f(x+1m)(=0,…,n-1 2 注:的确有局部却 即梯严需要2个初值和1来启动递推 但过程,这样的算法称为双步法 double-step 选、 method*/,而前面的三种算法都是单步法 / single-step method"1o a中点欧拉公式/m 中心差商近似导岁yG)s(x)rx)∠ 2h y(2)ay(y f(r,,y(,)) yH+1=y1+2hf(x,y)i=1,…,n-1 假设n=y(x11,男1=J(x),则可以导出R2=y(x1)-yn1=0(h) 即中点公式具有2阶精度
梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均 [ ( , ) ( , )] ( 0, ... , 1) 2 +1 = + f x y + f x +1 y +1 i = n − h y y i i i i i i 注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。 但注意到该公式是隐式公式，计算时不得不用到 迭代法，其迭代收敛性与欧拉公式相似。 ( ) ( ) 3 Ri = y xi+1 − yi+1 = O h 中点欧拉公式 /* midpoint formula */ 中心差商近似导数 h y x y x y x 2 ( ) ( ) ( ) 2 0 1 − x0 x1 x2 ( ) ( ) 2 ( , ( )) 2 0 1 x1 y x y x + h f x y yi+1 = yi−1 + 2h f (xi , yi ) i = 1, ... ,n−1 假设 ，则可以导出 即中点公式具有 2 阶精度。 ( ), ( ) i 1 i 1 i i y = y x y = y x − − ( ) ( ) 3 Ri = y xi+1 − yi+1 = O h 需要2个初值 y0和 y1来启动递推 过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步法 /* single-step method */
方法 显式欧拉 简单 精度低 隐式欧拉 稳定性最好精度低,计算量大 梯形公式 精度提高 计算量大 中点公式精度提高,显式 多一个初值, 可能影响精度 with OK,let’s edy y make it possible
方 法 显式欧拉 隐式欧拉 梯形公式 中点公式 简单 精度低 稳定性最好 精度低, 计算量大 精度提高 计算量大 精度提高, 显式 多一个初值, 可能影响精度 Can’t you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? OK, let’s Well, call me greedy… make it possible
改进欧拉法/ modified eulers method sepl:先用显式欧拉公式作预测,算出y+1=y1+hf(x,y) Sep2:再将代入隐式梯形公式的右边作校正,得到 y1=y1+-[f(x1,y)+f(x1+1,y+ #1=y;+[f(x,y)+∫(x1,y1+h∫(x1,y(i=0,…,n-1) 注:此法亦称为预测校正法/ predictor-corrector method+。 可以证明该算法具有2阶精度,同时可以看到它是个单 步递推格式,比隐式公式的迭代求解过程简单。后面将 看到,它的稳定性高于显式欧拉法
改进欧拉法 /* modified Euler’s method */ Step 1: 先用显式欧拉公式作预测，算出 ( , ) i 1 i i i y + = y + h f x y Step 2: 再将 yi+1 代入隐式梯形公式的右边作校正，得到 [ ( , ) ( , )] 2 +1 +1 +1 = + + i i i i i i f x y f x y h y y 注：此法亦称为预测-校正法 /* predictor-corrector method */。 可以证明该算法具有 2 阶精度，同时可以看到它是个单 步递推格式，比隐式公式的迭代求解过程简单。后面将 看到，它的稳定性高于显式欧拉法。 ( , ) ( , ( , )) ( 0, ... , 1) 2 +1 = + f x y + f x +1 y + h f x y i = n − h y y i i i i i i i i
§2龙格-库塔法/ Runge-Kutta Method G建立高精度的单步递推格式。 单步递推法的基本思想是从(x,y2)点出发,以某一斜 率沿直线达到(x1,y1)点。欧拉法及其各种变形所 能达到的最高精度为2阶 考察改进的欧拉法,可以将其改写为 斜率 y+1=y+列K1+K2 定取K1K2 的平均值吗? i K2 =f(i+h,vi+hK1) 步长一定是一个h吗?
§2 龙格 - 库塔法 /* Runge-Kutta Method */ 建立高精度的单步递推格式。 单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。 考察改进的欧拉法，可以将其改写为： ( , ) ( , ) 2 1 2 1 2 1 1 1 1 2 K f x h y hK K f x y y y h K K i i i i i i = + + = + = + + 斜率 一定取K1 K2 的平均值吗？ 步长一定是一个h 吗？
将改进欧拉法推广为: J+1 y2+h|K1+12K2 K f(xi, yu K2=f(x+ph, yi t pl y"(x)=,∫(x,y 首先希望能确定系数λ1、 精度,即在y;=y(x)的自 =f(x,y)+∫(x R;=y(x+1) d ∫(x,y)+∫1(x,y)f(x,y) Spl:将k2在(x1,y1)点作Tay K2=∫(x1+ph,y1+phK1) =f(x,y)+phf(x,y)+pMhKJf(x,y)+O2) y(x1)+py"(x;)+O(h2) Sep2:将K2代入第式,得到 Vam=V; +hay(x; )+a2ly'(x; )+phy(x, )+0(12)11 y;+(1+A2)hy(x;)+2ph2y"(x;)+O(h°)
首先希望能确定系数 1、2、p，使得到的算法格式有2阶 精度，即在 yi = y(xi ) 的前提假设下，使得 ( ) ( ) 3 Ri = y xi+1 − yi+1 = O h Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 ( , ) ( , ) ( , ) ( ) ( , ) 2 1 2 1 f x y phf x y phK f x y O h K f x ph y phK i i x i i y i i i i = + + + = + + ( ) ( ) ( ) 2 = y xi + phy xi + O h 将改进欧拉法推广为： ( , ) ( , ) [ ] 2 1 1 1 1 1 2 2 K f x ph y phK K f x y y y h K K i i i i i i = + + = + = + + ( , ) ( , ) ( , ) ( , ) ( , ) ( ) ( , ) f x y f x y f x y dx dy f x y f x y f x y dx d y x x y x y = + = + = Step 2: 将 K2 代入第1式，得到 ( ) ( ) ( ) ( ) ( ) [ ( ) ( ) ( )] 2 3 1 2 2 2 1 1 2 y hy x ph y x O h y y h y x y x phy x O h i i i i i i i i = + + + + + = + + + +
93:将与yxn)在点的泰勒展开作比较 y+1=y2+(x1+2)hy'(x;)+2ph2y"(x;)+O(h3) y(x+1)=y(x1)+hy(x1) n,(x)+O(h) 2 要求R1=y(x1)-yn1=O(h),则必须有: 1+A2=1,2P= 这里有3个未知 2 数,2个方程。 存在无穷多个解。所有满足上式的格式统称为2阶龙格-库 塔格式。注意到,p=1,1=就是叔进的欧拉法。 Q:为获得更高的精度,应该如何进一步推广?
Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 ( ) ( ) ( ) ( ) 2 3 yi+1 = yi + 1 + 2 hy xi + 2 ph y xi + O h ( ) ( ) 2 ( ) ( ) ( ) 3 2 1 y x O h h y xi+ = y xi + hy xi + i + 要求 ( ) ( ) ，则必须有： 3 Ri = y xi+1 − yi+1 = O h 2 1 1 , 1 + 2 = 2 p = 这里有 个未知 数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。 2 1 1, 注意到， p = 1 = 就是改进的欧拉法。 2 = Q: 为获得更高的精度，应该如何进一步推广？
p=y++2++km其中x(=1,…,m), K=f(i, vi) 2,…,m)和 K2=f(x+a2,y2+B21hk1) m;j=1,;-1)均为待定 系数,确定这些系数的 K=f(x+b+月1K+32h1)步骤与前面相似。 Km=f(x; +ah,y+Pm,hK+Bn hK,+.+hkrrd) >最常用为四级4阶经典龙格库塔法/ Classical Runge-Kutta Method * 1+=y+名(K1+2K2+2K3+K4) f(, yi) K2=∫(x+,y+2K1) K3=f(x2+,y2+K2) K = f(x; +h, y; +hk3
其中 i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i−1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。 ( , ... ) ...... ( , ) ( , ) ( , ) [ ... ] 1 1 2 2 1 1 3 3 31 1 32 2 2 2 21 1 1 1 1 1 2 2 − − + = + + + + + = + + + = + + = = + + + + m i m m m m m m i i i i i i i i m m K f x h y hK hK hK K f x h y hK hK K f x h y hK K f x y y y h K K K ➢ 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： ( , ) ( , ) ( , ) ( , ) ( 2 2 ) 4 3 3 2 2 2 2 2 2 1 1 1 6 1 2 3 4 K f x h y hK K f x y K K f x y K K f x y y y K K K K i i h i h i h i h i i i h i i = + + = + + = + + = + = + + + +