第九章常微分方程数值解 / Numerical Methods for Ordinary Differential Equations 6考虑一阶常微分方程的初值问题/ Initial-Value Problem*: d tf(x,y)x∈la,b y(a)=yo 只要∫(x,y)在{a,b×R上连续,且关于y满足 Lipschitz条 件,即存在与x,y无关的常数L使|(x,yn)-f(x,y2)SL|-y2 对任意定义在{a,b上的y(x)和y2(x)都成立,则上述IVP存 在唯一解。 要计算出解函数y(x)在一系列节点a=xx1<xn=b 处的近似值v≈y(x;)(=1,…,n) 节点间距h=x1-X(=…,n-1)为步长,通常采用等距节点, 即取h=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欧拉方注 亦称为欧拉折线法 向前去 Euler's polygonal arc method/ y(x)y/x(x)=y+b/(x,)运为y y+1=y2+hf(x1,y)(i=0,…,n-1) 定义在假设y=y(x),即第i步计算是精确的前提下,考 虑的截断误差R1=y(x1)-y称为局部截断误差 F local truncation error */ 定义若其然 阶精度。 R的主项OP"),则称该算法有p / leading term */ 欧拉法的局联在 R=y(x)h=1m)++号y(x)+O)+(川 ="y(x;)+O(h)欧拉法具有1阶精度
§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*/
8 1 Euler's Method >欧拉公式的改进: a隐式欧拉法/ implicit Euler method* 向后差商近似导数一y(x)yx)y(x yoyo thf(xu,y(D Vi+1=y; the Hey! Isnt the leading term of the loca truncation errorof Euler's method y(x)? 由于末为 Seems that we can make a good 到,故 称为隐式 use of it 式/ explicit*/ 欧拉公式。 一般先用显式计算一个初值,再迭代求解 隐式欧拉法的局部截断误差: R;=y(x11)-y=-2y(x)+O(h) 即隐式欧拉公式具有1阶精度
§1 Euler’s Method ➢ 欧拉公式的改进： 隐式欧拉法 /* 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
§1 Euler' s Method 梯形公式/ trapezoid formula一显、隐式两种算法的平均 J=2+,Uf(x,y)+f(x1m)(i=0,…,n-1 注:的确有局部都 即梯严需要2个初值y和y1来启动递推 但过程,这样的算法称为双步法/ double-step 迭 method,而前面的三种算法都是单步法 / single-step method*1o a中点欧拉公式/m 中心差商近似导岁y(x)=Xx)yx h →y(x2)≈y(yx2hf(x,y(x1) yH+1=y1+2hf(x,y)i=1,…,n-1 假设n=y(x1,y1=y(x),则可以导出R=yx-)-yn=0(h) 即中点公式具有2阶精度
§1 Euler’s Method 梯形公式 /* 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 */
&1 Euler's Method 方法 显式欧拉 简单 精度低 隐式欧拉 稳定性最好精度低,计算量大 梯形公式 精度提高计算量大 中点公式精度提高,显式 多一个初值, 可能影响精度 with OK,let’s ledy make it possible
方 法 §1 Euler’s Method 显式欧拉 隐式欧拉 梯形公式 中点公式 简单 精度低 稳定性最好 精度低, 计算量大 精度提高 计算量大 精度提高, 显式 多一个初值, 可能影响精度 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
§1 Euler' s Method a改进欧拉法/ modified eulers method Se1:先用显式欧拉公式作预测,算出y=y1+hf(x1,y;) Se2:再将代入隐式梯形公式的右边作校正,得到 y1=y1+[∫(x;,y)+f(x+1,y 2 #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 §1 Euler’s Method
§2龙格-库塔法/ Runge-Kutta Method 建立高精度的单步递推格式。 单步递推法的基本思想是从(x,y)点出发,以某一斜 率沿直线达到(x1,y)点。欧拉法及其各种变形所 能达到的最高精度为2阶。 6考察改进的欧拉法,可以将其改写为: 斜率 y+1=y+列K1+K2 一定取K1K2 的平均值吗? f(x,y K2 =f(i+h,vi+hKi 步长一定是一个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 吗？
§2 Runge- Kutta method 将改进欧拉法推广为: y+=y1+hK1+2K2 K f(xi, yu K =f(r +ph, y;+plr y(x)=,f(x, y) 首先希望能确定系数λ1 精度,即在y;=y(x)的自 =f(x,y)+∫(x R=y(x+1) d ∫(x,y)+∫1(x,y)f(x,y) sepl:将k2在(x1,y2)点作Tay K2=f(x+ph, y; phU) f(xi,y)+phf,(x;, yi+phK,(i,yi)+O(h) -y(x,)+phy(x1)+O() Sp2:将K2代入第式,得到 Va=D2 +hay'(x;)+aly'(x )+ phy(x)+O(2) y1+(1+2)hy(x;)+a2ph2y"(x;)+O(h3)
§2 Runge-Kutta Method 首先希望能确定系数 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 = + + + + + = + + + +
§2 Runge- Kutta method Sep3:将与y(x1)在x点的泰勒展开作比较 y+1=y+(1+A2)hy'(x;)+2ph2y"(x)+O(h) y 0=y(xi)+hy(i )+y(x)+o(h 2 要求R1=y(x-1)-y1=O(h),则必须有: 1+A2=1,2P= 这里有3个未知 2 数,2个方程。 存在无穷多个解。所有满足上式的格式统称为2阶龙格-库 塔格式。注意到,p=,41=就是改进的欧拉法。 Q:为获得更高的精度,应该如何进一步推广?
§2 Runge-Kutta Method 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: 为获得更高的精度，应该如何进一步推广？
§2 Runge- Kutta method h=+川1+K2+,+nkm其中λ(i=1,…,m),a1( K=f(,vi) 2,…,m)和(i=2,…, K2=f(x+2b,y2+B21hk1) m;j=1,;-1)均为待定 K3=八(x+a4+3+2h)系数,确定这些系数的 步骤与前面相似。 K=f(x; +am h,y+Bm, hK+B hK,.+BhK) 最常用为四级4阶经典龙格库塔法/ Classical Runge-Kutta Method * y+1=y1+b(K1+2K2+2K3+K4) =f(x1,y) K2=∫(x1+,y+K1) K3=f(x1+,y1+2K2) f(+h, y; +hk3)
其中 i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i−1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。 §2 Runge-Kutta Method ( , ... ) ...... ( , ) ( , ) ( , ) [ ... ] 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 = + + = + + = + + = + = + + + +