常微分方程的求解.pdf

18— 1 常微分方程数值解法 2§ 1 引言§ 2 Euler方法§ 3 Runge－ Kutta方法§ 4 单步法的收敛性与稳定性§ 5 线性多步法§ 6 方程组与高阶方程的情况§ 7 边值问题的数值解法3§1 引言 微分方程 ：关于一个未知函数的方程，方程中含有未知函数的（偏）导数，以及自变量等，其中关于未知函数导数的最高次数称为微分方程的阶数 . 例如：0)()(')()('' =++− xcyxbyxaxy4实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在理论研究与工程实际上应用很广泛. 很多问题的数学模型都可以归结为常微分方程. 很多偏微分方程问题，也可以化为常微分方程问题来近似求解. 微分方程的 应用情况 5对于一个常微分方程： '(,) , [,]dyyfxyxabdx== ∈为了使解存在，一般要对函数 f 施加限制条件，例如要求 f 对 y 满足 Lipschitz 条件： 1212(, ) (, )f xy fxy Ly y−≤−6同时， 一个有解的微分方程通常会有无穷多个解例如 cos( ) sin( ) ,dyx yxaaRdx= ⇒= +∀∈为了使解唯一，需要加入一个限定条件. 通常会在端点出给出，如下面的 初值问题 ： 0(, ), [,]()dyf xy x abdxya y⎧=∈⎪⎨⎪=⎩7常微分方程的解是一个函数，但是，只有极少数特殊的方程才能求解出来，绝大多数是不可解的.并且计算机没有办法对函数进行运算. 一般考虑其近似解法，一种是近似解析法，如逼近法、级数解法等，另一种是本章介绍的 数值解法 . 8§2 Euler 方法 92－ 1 Euler 公式 对常微分方程初值问题： ⎩⎨⎧==00')(),(yxyyxfy数值求解的关键在于消除其中的导数项——称为 离散化 . 利用差商近似逼近微分是离散化的一个基本途径 . 10现在假设求解节点为),,1,0( miihaxi"=+=，其中mabh−=为步长，这些节点相应的函数值为)(,),(1 mxyxy ". 在点 nx处，已知 ))(,()('nnnxyxfxy =用 nx的向前差商nnnnxxxyxy−−++11)()(近似代替)('nxy，如§ 1，则得到所谓的 Euler 公式 1(,)nn nnyyhfxy+= + —— 单步、显式格式 11Euler 公式的局部截断误差： 假设)(nnxyy =情况下， 11)(++−nnyxy称为 局部截断误差 . '''2311''23()( ) () () ()2() (,()(()))2nnn n nnnnnyxyx y yx hy x h Ohyx hyxfxyx hOh++−= + + +−− =+故有)(2)(''211 nnnxyhyxy ≈−++ . 122－ 2 后退的 Euler 公式 同样对常微分方程初值问题，在 1+nx点，已知))(,()(111'+++=nnnxyxfxy，如果用向后差商hxyxynn)()(1−+代替)(1'+nxy， 则得到 后退的 Euler 公式： 111(,)nn nnyyhfxy+ ++= +—— 单步、隐式格式 13相对于以上可以直接计算 1+ny的 Euler 公式（显式） ，上式是隐式公式 . 一般来讲，显式容易计算，而隐式具有更好的稳定性 . 求解上述公式，通常使用 迭代法 ：对于给定的初值)0(1+ny，计算 (1) ()111(, (0,1,)kknnnnyyfxyk++++=+ = "， 如果)(1limknky+∞→收敛，则其极限必满足上述后退Euler 公式 . 14局部截断误差 ： 假设)(nnxyy =，则 ),()(111 ++++=nnnnyxhfxyy. 由于 )]()[,())(,(),(1111111 +++++++−+=nnnynnnnxyyxfxyxfyxf η且 ''' 211 1(,() () () () ()nn n n nf xyx yx yx hyxOh++ +==++15则有 '2' 31111(,)[ ()]() () () ()nyn n n n n nyhfx yyx yxhyxhyxOhη++++=−+++将此式减去式 2'''31( ) () () () ()2nnn nhyx yx hy x y x Oh+=+ + +可得， 2'' 311 1 11() (,)[() ] ()()2nnyn nn nhy xyhfx yxy yxOhη++ + ++−= −− +16考虑到21111(,)()1(,)ynynhf xOhhf xηη++=+ +−，则有 22'' 3 ''11( ) () () ()nn n nhhy x y yx Oh yx++−=− + ≈−172－ 3 梯形公式 由于上述两个公式的局部截断误差绝对值相等，符号相反，故求其算术平均得到 梯形公式 ： 111[( , ) ( , )]2nn nn nnhyy fxyfxy+++=+ +—— 单步、隐式格式 18梯形法同样是隐式公式，可用下列迭代公式求解：(0)1(1) ()111(, )[( , ) ( , )]2nn nnkknn nn nnyyhfxyhyyfxyfxy+++++⎧ =+⎪⎨=+ +⎪⎩局部截断误差：类似于后退 Euler，可计算出)(12)('''311 nnnxyhyxy −≈−++192－ 4 改进的 Euler 公式 上述用迭代法求解梯形公式虽然提高了精度，但计算量也很大 . 实际上常采用的方法是， 用 Euler公式求得初始值（预测） ，然后迭代法仅施行一次（校正）—— 改进的 Euler 公式： 1111(,)[( , ) ( , )]2nnnnnn nn nnyyfxyhyy fxyfxy++++⎧=+⎪⎨=+ +⎪⎩20估计上式中第二式当 1+ny为准确值时的局部截断误差： ''11 1 13(3)() ()(() [() ()])2()12nn n n n nnhyx y yx yx y x y xhyx++ + +−= − + +≈−212－ 5 Euler 两步公式 如果用中心差商hxyxynn2)()(11 −+−代替)('nxy， 则得 Euler 两步公式 112(,)nn nnyy hfxy+−=+—— 两步、显式格式22假设 1−ny及 ny均为准确值，利用 Taylor 展式容易计算 Euler 两步公式的 局部截断误差 为：11 1 13(3)() ()(()2(,())()3nn n n nnnyx y yx yx hf x yxhyx++ + −−= − +≈23此式与梯形公式相结合，得到如下的 预测－校正公式： 111112(,)[( , ) ( , )]2nnnnnn nn nnyy hfxyhyy fxyfxy−++++⎧=+⎪⎨=+ +⎪⎩假设第一式中的 1−ny及 ny，以及第二式中的 ny及1+ny均是准确值，则有， 2441)()(1111−≈−−++++nnnnyxyyxy从而可得以下的事后估计式， 111111 114() ( )51() ( )5nnnnnn nny x yyyyx y y y+ +++++ ++⎧−≈− −⎪⎪⎨⎪−≈ −⎪⎩25可以期望，以上式估计的误差作为计算结果的补偿，可以提高计算精度 . 以 np及 nc分别表示第 n 步的预测值和校正值， 则有以下的“ 预测－改进－校正－改进 ”方案（其中在 1+np与 1+nc尚未计算出来的前提下，以nncp −代替 11 ++−nncp： 26预测：'112nnnhyyp +=−+ 预测的改进：)(5411 nnnncppm −−=++ 计算：),(11'1 +++=nnnmxfm校正：)(2'1'1 ++++=nnnnmyhyc校正的改进： )(511111 ++++−+=nnnncpcy计算：),(11'1 +++=nnnyxfy27例 用 Euler 方法求解初值问题 2'[0,0.6](0) 1yyxyxy⎧ =− −∈⎨=⎩取0.2h=，要求保留六位小数 . 解： Euler 迭代格式为 2210.2( ) 0.8 0.2kk kkk k kkyy yxy y xy+=+ −− = −因此 282100(0.2) 0.8 0.2 0.8yyyxy≈= − =2211(0.4) 0.8 0.2 0.6144yyyxy≈= − =232(0.6) 0.8 0.2 0.461321yyyxy≈= − =29例 用改进的 Euler 方法求解初值问题 2'sin0[0,0.6](0) 1yyy xxy⎧ ++ =∈⎨=⎩取 0.2h= ，求(0.2), (0.4)yy的近似值，要求保留六位小数 . 解：改进的 Euler 格式为 212211110.2( sin )0.2(sin sin)2kkkkkkk kkk kkkyy yyxy yyyxyyx+++++⎧ =+ −−⎪⎨=+ −− − −⎪⎩30即， 2221 10.82 0.08 sin 0.1(0.8 0.2 sin ) sinkkkk kkkkyyyx yyxx+ +=− − −则有 1(0.2) 0.807285yy≈ =， 2(0.4) 0.636650yy≈=31§3 Runge－ Kutta 方法 Def.1 如果一种方法的局部截断误差为)(1+phO，则称该方法具有 p 阶精度 . 323－ 2 Runge— Kutta 方法的基本思想 上述的 Taylor 级数法虽然可得到较高精度的近似公式，但计算导数比较麻烦 . 这里介绍不用计算导数的方法 . ))(,()()()('1hxyhxfhxyhxyxynnnnnθθθ ++=+=−+—— 平均斜率 . 33如果粗略地以),(nnyxf作为平均斜率，则得Euler 公式； 如果以221KK +作为平均斜率，其中),(1 nnyxfK =，),(112hKyxfKnn+=+ ，则得改进的 Euler 公式 . 343－ 3 二阶的 Runge－ Kutta 方法 对点 nx和)10( ≤<+=+pphxxnpn ，用这两点斜率的线性组合近似代替平均斜率，则得计算公式：1112121()(, )(, )nnnnnp nyyhKKKfxyKfxyphKλλ++⎧ =+ +⎪=⎨⎪=+⎩35现确定系数p,,21λλ，使得公式具有二阶精度 . 因为，取 ny为()nyx，则 '1(, ) (,() '()nn n n n nKfxy fxyx yx y== =再把 2K在),(nnyx处展开，有 36'21(, )( , )np n n n nK f x y phK f x ph y phy+=+=++代入可得， '2'31122() ()nn n ny yhyphyOhλλ λ+=++ + +'2(, ) (, ) (, ) ()nn xnn ynn nf xy fxy ph fxy phy Oh=+⋅+⋅+'2(')(,)()nxynnyphffyxyOh=+⋅+⋅ +'''2()nnyphyOh=+⋅+37相比较二阶 Taylor 展开''2'12nnnnyhhyyy ++=+ ，有，⎪⎩⎪⎨⎧==+211221pλλλ满足此条件的公式称为 二阶 Runge－ Kutta 公式 .38可以验证改进的 Euler 公式属于二阶 Runge－Kutta 公式 . 下列变形的 Euler 公式也是二阶Runge－ Kutta 公式： 12121(, )(, )22nnnnnnyyhKKfxyhhK fx y K+⎧⎪=+⎪=⎨⎪⎪ =++⎩393－ 4 三阶 Runge－ Kutta 公式。