NA005a微分方程求解

华中科技大学 李红 第五章 常微分方程数值解 /* Numerical Methods for Ordinary Differential Equations */ �考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: ⎪⎩ ⎪⎨ ⎧ = ∈= 0)( ],[),( yay baxyxf dx dy 只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述IVP存 在唯一解。 |||),(),(| 2121 yyLyxfyxf −≤− 本章的任务：计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值。),...,1()( Nnxyy nn =≈ 华中科技大学 李红 建立常微分方程数值方法的基本思想 微分方程数值解法，其实是求出方程的解 在一系 列离散点上的近似值。则微分方程数值解的基本思想 是：求解区间和方程离散化。 )(xy 将求解区间 离散化,是在 上插入一系列的分 点 ,使 bxxxxxa Nnn =<<<<<<= + ...... 110 [ ]ba, [ ]ba , { }kx ， 求解区间离散化¾ 记 称为步长，一般取hn = h 称为等步长节点 。 )1,...,1,0(1 −=−= + Nnxxh nnn nhxxn += 0 ),...,2,1,0( Nn = N abh −=(常数),节点为 华中科技大学 李红 ∫=− m n x xnm dxxyxfxyxy ))(,()()( ))(( 00 yxy = 然后将上式右端采用第四章介绍的数值积分离散化， 从而获得原初值问题的一个离散差分格式。 将微分方程离散化¾ 将微分方程离散化，通常有下述方法： （1）差商逼近法 即是用适当的差商逼近导数值。 （2）数值积分法 基本思想是先将问题转化为积分方程 （3）Taylor展开法 见后面的叙述。 华中科技大学 李红 §1 欧拉方法 /* Euler’s Method */ ¾欧拉公式： x0 x1 向前差商近似导数 h xyxyxy )()()( 010 −≈′ ),()()()( 000001 yxfhyxyhxyxy +=′+≈ 1y记为 )1,...,0(),(1 −=+=+ Nnyxfhyy nnnn ⎪⎩ ⎪⎨ ⎧ = −=′ 1)0( 2 y y xyy )10( << x 例. 求初值问题 解：本题的Euler公式的具体形式为 )( n n nnn y xyhyy 21 −+=+ §1 Euler’s Method 华中科技大学 李红 1.73211.78481.01.41421.43510.5 1.67331.71780.91.24161.35820.4 1.61251.64980.81.26491.27740.3 1.54921.58030.71.18321.19180.2 1.48321.50900.61.09541.10000.1 nx nxny ny)( nxy )( nxy 取步长h=0.1。我们将计算结果与其解析解的精确值一同 列在下中，其中 是节点， 是节点上的近似值， 是精确值，结果见下表： nynx )( nxy 如果说表格仍不够直观的话，我们用Matlab做出积分曲 线与近似值的图，如下： §1 Euler’s Method 华中科技大学 李红 从图中可以看出，灰色连续的曲线就是初值问题的解析解 的曲线。而红色的星点便是数值解。在图上似乎 数值解与曲线的偏差不是很大，但不要忘记这只是在0到1范围 内的。通过后面用其他方法解本题，大家便会发现Euler方法误 差其实是很大的。 xy 21 += Matlab 作图显示 §1 Euler’s Method 华中科技大学 李红 §1 Euler’s Method # 隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 h xyxyxy )()()( 011 −≈′ x0 x1 ))(,()( 1101 xyxfhyxy +≈ ),,( =+y 1 (+ fh 1 1...0, 1+ Nnx n −= yy nnn )+ 由于未知数 yn+1同时出现在等式的两边， 不能直接得到，故称为隐式 /* implicit */ 欧拉公式，而前者称 为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。即 ),()0( 1 nnnn yxhfyy +=+ ),( )( 11 )1( 1 k nnn k n yxhfyy ++ + + += L,2,1,0=k 如果迭代过程收敛，则某步后 就可以作为 ，从而进行 下一步的计算。 )( 1 k ny + 1+ny 华中科技大学 李红 梯形方法的平均化思想 可以借助几何直观说明，同 Euler方法的图，见右： §1 Euler’s Method# 梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均 )1,...,0()],(),([ 2 111 −=++= +++ Nnyxfyxfhyy nnnnnn # 两步欧拉公式 /* midpoint formula */ 中心差商近似导数 h xyxyxy 2 )()()( 021 −≈′ x0 x2x1 ))(,(2)()( 1102 xyxfhxyxy +≈ 1,...,1),(211 −=+= −+ Nnyxfhyy nnnn nP A B nx 1+nx Q 1+nP 华中科技大学 李红 # 改进欧拉法 /* modified Euler’s method */ Step 1: 先用显式欧拉公式作预测，算出 ),(1 nnnn yxfhyy +=+ Step 2: 再将 代入隐式梯形公式的右边作校正，得到1+ny )],(),([ 2 111 +++ ++= nnnnnn yxfyxfhyy 注：此法亦称为预测-校正法 /* predictor-corrector method */。 可以证明该算法具有 2 阶精度，同时可以看到它是个单 步递推格式，比隐式公式的迭代求解过程简单。后面将 看到，它的稳定性高于显式欧拉法。 ( )[ ] )1,...,0(),(,),( 2 11 −=+++= ++ Nnyxfhyxfyxfhyy nnnnnnnn §1 Euler’s Method )( 2 1 ),(,),( 1 1 cpn pnncnnnp yyy yxhfyyyxhfyy += +=+= + + 华中科技大学 李红 解：改进的Euler公式为 ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ += −+= −+= + + )( )( )( cpn p n pnc n n nnp yyy y xyhyy y xyhyy 2 1 2 2 1 1 例 作为比较，我们仍用Euler法中的那个例子。 ⎪⎩ ⎪⎨ ⎧ = −=′ 1)0( 2 y y xyy )10( << x 1.73211.73791.0 1.67331.67820.9 1.61251.61530.8 1.54921.55250.7 1.48321.48600.6 1.41421.41640.5 1.34161.34340.4 1.26491.26620.3 1.18321.18410.2 1.09541.09590.1 nx ny )( nxy 我们仍取步长h=0.1，计算结果 见右表： §1 Euler’s Method 华中科技大学 李红 Matlab作图显示 从表和图可以看出，改进的Euler法的精度提高了 不少。 §1 Euler’s Method 华中科技大学 李红 初值问题的单步法可用一般形式表示为: ¾ 局部截断误差和方法的阶 ),,,( 11 hyyxhyy nnnnn ++ += ϕ 称为增量函数ϕ ),(),,,( 111 +++ = nnnnn yxfhyyxϕ 隐式欧拉法有例如对欧拉法有 ，),(),,( nnnn yxfhyx =ϕ 含有 时，方法是隐式的， 不含有 时，方法是显式的。 1+nyϕϕ 1+ny 从 开始计算，如果考虑每一步产生的误差，直到 则有误差 ，称为该方法在 的整体截断误 差，和求得整体截断误差是复杂的。为此，我们仅考 虑从 到 的局部情况，并假设 之前的计算没有误 差，即 ,下面给出单步法的局部截断误差概念。 nnn yxye −= )( 0x nx nx nx 1+nx nx )( nn xyy = §1 Euler’s Method 华中科技大学 李红 即 在假设 yn = y(xn)，即第n 步计算是精确的前提下， 考虑的截断误差 Tn+1= y(xn+1) − yn+1 称为局部截断误差 /* local truncation error */。 定义 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。 )欧拉法的局部截断误差： )],([)]()()()([)( 32111 2 nnnn h nnnnn yxhfyhOxyxyhxyyxyT +−+′′+′+=−= +++ )()( 32 2 hOxy nh +′′= 欧拉法具有 1 阶精度。 为显式单步法的局部截断误差。 )),(,()()( 11 hxyxhxyxyT nnnnn ϕ−−= ++ 定义 设 是初值问题的准确解，称)(xy 即两步欧拉公式具有 2 阶精度。)()( 3111 hOyxyT nnn =−= +++ )两步欧拉法的局部截断误差： §1 Euler’s Method 华中科技大学 李红 )()( 2 )]()()([)()( 2 )( ))(,()()( 3 2 23 2 1111 hxyh hxyhxyhhxyhxyh xyxhfxyxyT n nnnn nnnnn Ο+′′−= Ο+′′+′−Ο+′′+′= −−= ++++ 例 求隐式欧拉格式 的局部截断误差。),( 111 +++ += nnnn yxhfyy 具有 1 阶精度。 例 求梯形格式 的局部 截断误差。 )],(),([ 2 111 +++ ++= nnnnnn yxfyxfhyy 以上定义对隐式单步法也是适用的。 )()( 12 )()]( 2 )()()([ 2 )( !3 )( 2 )( )]()([ 2 )()( 4 3 4 2 32 111 hxyh hxyhxyhxyxyh xyhxyhxyh xyxyhxyxyT n nnnn nnn nnnnn Ο+′′′−= Ο+′′′+′′+′+′− ′′′+′′+′= ′+′−−= +++ 具有 2 阶精度。 §1 Euler’s Method 华中科技大学 李红 显然， , 是在点 , 处的斜率。以上公式用到 的是两个点的斜率的加权平均，它为构造算法提供了新 的途径。Runge-Kutta方法就是这种思想的体现和发展。 1k 2k nx 1+nx §2 龙格 - 库塔法 /* Runge-Kutta Method */ 建立高精度的单步递推格式。 单步递推法的基本思想是从 ( xn , yn ) 点出发，以某一 斜率沿直线达到 ( xn+1 , yn+1 ) 点。欧拉法及其各种变形 所能达到的最高精度为2阶。 �考察改进的欧拉法，可以将其改写为： ),( ),( 2 1 2 1 12 1 211 hKyhxfK yxfK KKhyy nn nn nn ++= = ⎥⎦ ⎤⎢⎣ ⎡ ++=+ 斜率 一定取K1 K2 的平均值吗？ 步长一定是一个h 吗？ §2 Runge-Kutta Method ))(,()()( ),(),()()( 1 1 1 ξξ ξξ yhfxyxy xxy h xyxy nn nn nn += ∈′=− + + + 华中科技大学 李红 §2 Runge-Kutta Method 首先希望能确定系数 λ1、λ2、p，使得到的算法格式有2阶 精度，即在 的前提假设下，使得 )()( 3111 hOyxyT nnn =−= +++ )( nn xyy = Step 1: 将 K2 在 ( xn , yn ) 点作 Taylor 展开 )(),(),(),( ),( 2 1 12 hOyxfphKyxphfyxf phKyphxfK nnynnxnn nn +++= ++= )()()( 2hOxyphxy nn +′′+′= Step 2: 将 K2 代入第1式，得到{ } )()()()( )]()()([)( 32 221 2 211 hOxyphxyhy hOxyphxyxyhyy nnn nnnnn +′′+′++= +′′+′+′+=+ λλλ λλ 将改进欧拉法推广为： ),( ),( ][ 12 1 22111 phKyphxfK yxfK KKhyy nn nn nn ++= = ++=+ λλ 121 =+ λλ 10 ≤< p 华中科技大学 李红 §2 Runge-Kutta Method Step 3: 将 yn+1 与 y( xn+1 ) 在 xn点的泰勒展开作比较 )()()()( 322211 hOxyphxyhyy nnnn +′′+′++=+ λλλ )()( 2 )()()( 3 2 1 hOxy hxyhxyxy nnnn +′′+′+=+ ，则必须 有： )()( 3111 hOyxyT nnn =−= +++ 2 1,1 221 ==+ pλλλ 这里有 个未知数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。 2 1,1 21 === λλp注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？ 华中科技大学 李红 ),(1 nn yxfk = .............. ∑ = + += L i iinn khyy 1 1 λ ))(,( 232131333 kakahcyhcxfk nn +++= 其中 , , 均为待定系数，确定这些系数 的步骤与前面相似。 1 1 =∑ = L i iλ ∑− = = 1 1 1 i j ija1≤ic ),( 1 1 ∑− = ++= i j jijinini kahcyhcxfk Li ,,3,2 L= ),( 1222 hkcyhcxfk nn ++= Runge-Kutta方法的一般形式 §2 Runge-Kutta Method 华中科技大学 李红 §2 Runge-Kutta Method ¾最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK KyxfK yxfK KKKKyy nn h n h n h n h n nn h nn ++= ++= ++= = ++++=+ ⎪⎩ ⎪⎨ ⎧ = −=′ 1)0( 2 y y xyy )10( << x 我们仍用前面的例子来看看四阶Runge-Kutta方法的效果。 华中科技大学 李红 解: 对于本题，经典的四阶Runge-Kutta方法具有以下形式： ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ + +−+= + +−+= + +−+= −= ++++=+ 3 34 2 23 1 12 1 43211 )(2 2 2 2 2 2 2 2 )22( 6 hky hxhkyk khy hxkhyk khy hxkhyk y xyk kkkkhyy n n n n n n n n n n n n nn §2 Runge-Kutta Method 1.73211.73211.0 1.61251.61250.8 1.48321.48330.6 1.34161.34170.4 1.18321.18320.2 nx ny )( nxy 这里，我们取步长 h=0.2，下面是计算结果 （符号的意义同前） 华中科技大学 李红 从上图可以看出，每一个数值解都“准确”的落在 了真实解的曲线上。与改进的Euler法所做出来的图比 较，似乎精确性没有很明显的提高，但是不要忘了在 用Runge-Kutta方法时，我们取的步长是h=0.2。实际 上，Runge-Kutta方法的精确性是要高很多。 §2 Runge-Kutta Method Matlab作图 华中科技大学 李红 §2 Runge-Kutta Method 注： )龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的 值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系： 753 可达到的最高精度 642每步须算Ki 的个数 )( 2hO )( 3hO )( 4hO )( 5hO )( 6hO)( 4hO )( 2−nhO 8≥n ) 由于龙格-库塔法的导出基于泰勒展开，故精度主要受 解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。 HW: p.147-148 #2, # 3, # 4 华中科技大学 李红 §3 收敛性与稳定性 /* Convergency and Stability */ ¾ 收敛性 /* Convergency */ 定义 若某算法对于任意固定的 x = xn= x0 + n h，当 h→0 ( 同时 n→ ∞) 时有 yn→ y( xn )，则称该算法是收敛的。 例：就初值问题 考察欧拉显式格式的收敛性。 ⎩⎨ ⎧ = =′ 0)0( yy yy λ 解：该问题的精确解为 xeyxy λ0)( = 欧拉公式为 nnnn yhyhyy )1(1 λλ +=+=+ 0)1( yhy nn λ+= 对任意固定的 x = xn = n h ，有 nn xhhx n hyhyy λλλλ ])1[()1( /10/0 +=+= )(0 n x xyey n =λ→ 9 §3 Convergency and Stability 华中科技大学 李红 注：1、判断显式单步格式的收敛性，归结为验证增量函数 是否满足 Lipschitz 条件。 2、单步格式的整体截断误差由初值误差及局部截断误差决定，整 体截断误差比局部截断误差的阶数低一阶。 3、要构造高精度的计算方法，只需设法提高局部截断误差的阶即 可。 ϕ （2）微分方程的初值是精确的. 。)()( pnnn hyxye Ο=−= 对于一个p 阶的显式单步法，若满足如下条件定理 （1）增量函数 关于y 满足Lipschitz条件，即存在常数 0>ϕL ，使 ϕ RyyyyLhyxhyx ∈∀−≤− ,,),,(),,( ϕϕϕ 成立； §3 Convergency and Stability 更详细的请参看教材p132 则该方法收敛，其整体截断误差为: 华中科技大学 李红 §3 Convergency and Stability ¾稳定性 /* Stability */ 例：考察初值问题 在区间[0, 0.5]上的解。 分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 ⎩⎨ ⎧ = −=′ 1)0( )(30)( y xyxy 0.0 0.1 0.2 0.3 0.4 0.5 精确解改进欧拉法欧拉隐式欧拉显式节点 xi xey 30−= 1.0000 −2.0000 4.0000 −8.0000 1.6000×101 −3.2000×101 1.0000 2.5000×10−1 6.2500×10−2 1.5625×10−2 3.9063×10−3 9.7656×10−4 1.0000 2.5000 6.2500 1.5626×101 3.9063×101 9.7656×101 1.0000 4.9787×10−2 2.4788×10−3 1.2341×10−4 6.1442×10−6 3.0590×10−7 华中科技大学 李红 §3 Convergency and Stability 定义 若某算法在计算过程中任一步产生的误差在以后的计 算中都逐步衰减，则称该算法是绝对稳定的 /*absolutely stable */。 一般分析时为简单起见，只考虑试验方程 /* test equation */ 0Re( <=′ )λλ yy 常数，可以 是复数 当步长取为 h 时，将某算法应用于上式，并假设只在初值 产生误差 ，则若此误差以后逐步衰减，就称该 算法相对于 绝对稳定， 的全体构成绝对稳定区 域。我们称算法A 比算法B 稳定，就是指 A 的绝对稳定区域 比 B 的大。 000 yy −=ε h λ h= h 华中科技大学 李红 §3 Convergency and Stability 例：考察显式欧拉法 011 )1( yhyhyy iiii ++ +=+= λ 000 yy −=ε 011 )1( yhy ii ++ += 0 1 111 )1( εε ++++ +=−= iiii hyy 由此可见，要保证初始误差ε0 以后逐步衰减， 必须满足： hh λ= 1|1| <+ h 0-1-2 Re Img 例：考察隐式欧拉法 11 ++ += iii yhyy λ ii yh y ⎟⎠ ⎞⎜⎝ ⎛ −=+ 1 1 1 0 1 1 1 1 εε + + ⎟⎠ ⎞⎜⎝ ⎛ −= i i h 可见绝对稳定区域为： 1|1| >− h 210 Re Img 注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法 的好。 华中科技大学 李红 §3 Convergency and Stability 例：隐式龙格-库塔法 ⎪⎩ ⎪⎨ ⎧ = ++++= +++=+ ),...,1( )...,( ]...[ 11 111 mj hKhKyhxfK KKhyy mmjjijij mmii ββα λλ 而显式 1~ 4 阶方法的绝对稳定 区域为 ⎪⎩ ⎪⎨ ⎧ ++= +=+ ) 2 , 2 ( 11 11 KhyhxfK hKyy ii ii 其中2阶方法 的绝对稳定区域为 0 Re Img k=1 k=2 k=3 k=4 -1-2-3 - - - 1 2 3 Re Img 无条件稳定 HW: p.148 #6,8 建立常微分方程数值方法的基本思想 Matlab作图显示