第九章+常微分方程初值问题数值解法

第九章常微分方程初值问题

的数值解法在自然科学的工程技术的许多领域中，常会遇到常微分方程初值问题，但这种问题大多数情况下不存在初等形式的解析解，只能用近似方法来求解。近似解法主要有两类：一类叫做近似解析方法，它能给出解得近似表达式，例如熟知的级数解法和逐次逼近法等；另一类近似解法称为数值解法，它可以给出解在一些离散点上的近似值。数值方法便于电子计算机求解微分方程。本章讨论初值问题常用的数值解法，并介绍有关的基本理论。假设给定一阶常微分方程的初值问题：dy(9.1)(9.2)<dx=f(x，y)，a<x<b、(9.1)(9.2)其中f为x,y的已知函数，a为给定的初始值。由微分方程的理论可知，如果f{x,y)在区域a<x<b即存(9.2)则是-c»vyv+«内连续，且关于y满足李普希兹(Lipschitz)条件，在常数即存(9.2)则是If(x,y)-f(x,y)<LIy一yl对所有的a<x<b及任何y，y均成立，则初值问题(9.1)，有连续可微的解y(x)存在且唯一。所谓初值问题的数值解，问题的解y(x)在一系列点a=x0vxivx2v...vx1vx2=b处的值y(x)的近似值y(n=0,1,・・・,N)。这里相邻两个节点之间的距离h=x\-x通常称为步长，通常将步长h取为常数h。§9.1欧拉方法与改进的欧拉方法欧拉(Euler)方法是最简单的数值方法，由于它的精确程度较差，已不常用于实际计算。但构造这个方法的基本原理，对于构造一般的数值方法具有普遍意义，因此首先对它进行讨论。1欧拉方法的构造将方程(9.1)中点x处的导数 y(x)用差商近似地表示为TOC\o"1-5"\h\zy(x)^-y(x)y(x)a -—八- h即在该点有近似等式y(、4y(七)af(x,y(x))h - -用近似值y代替y(x)，由上式可以导出其近似值满足的差分方程：y1=y+hf(x,y(x)), n=0,1,...,N-1 (9.3)(9.3)式称为欧拉方法的计算公式或称为欧拉公式。当初始值y=a给定时，利用欧拉公式就可以逐次计算出初值问题的数值0解*，yLyN。2后退的欧拉公式如果在点x处用向后差商近似(9.1)式中的导数-+1y(x)一y(x)y(x)a —+1 n_- h即有 y(\+1)ay(xp+hf(x-+『y(x〃,「)再用近似值y.代替y(x.)0=-,-+1)，则可导出近似值满足的差分方程

yn+1=yn+hf(x^+1，七+1) (9.6)(9.6)式称为后退的欧拉公式。这是一个两端含有未知数y的n+1方程，这样的方法称为隐式方法。相应地，称右端不隐含y的n+1欧拉公式(9.3)是一个显式方法。隐式方法求y时需解方程，n+1当f(X,y)关于y为非线性时，就需要解非线性方程，通常用迭代法求解。可以证明，隐式的后退欧拉公式(9.6)的局部截断误差为(9.7)R=-兽y〃(x)+O伽)n,h 2n(9.7)事实上，设y=y(x)时，由(9.6)式有y=y+hf(x,y)(9.8)n+1n n+1n+1(9.8)=y(x)+hf(x,y)n n+1n+1由中值定理得(9.9)f(xn+1，yn+1)(9.9)TOC\o"1-5"\h\z=f(x,y(x))+f'(x疝)[y-y(x)]n+1 n+1 yn+1 n+1 n+1得(9.10)将(9.9)式代入(9.8)式，并利用微分方程得(9.10)y =y(x)+hyf(x)n+1 n n+1+hf'(x,n)[y-y(x)]yn+1 n+1 n+1由微分方程解的泰勒展式(9.4)与(9.10)式相减，得\o"CurrentDocument"y(x)-yn+1 n+1h2=h[y(x)-y(x)]+ y(x)-hf(x户)n n+1 2nyn+1.[y-y(x)]+o(h3)n+1 n+1以yXX)=yxx)+hyf,(x)+O(h3)代入上式，得TOC\o"1-5"\h\zn+1 n ny(Xn+1)-yn+1h2 一=-y(X)+hf(X如)[y(X)-y]+O((h3)2 nyn+1 n+1 n+1移项整理得y(x)-yn+1 n+1= 1 [-性y〃(X)+o(h3)]1-hf(X,n) 2nyn+1注意当h充分小，使时;|<1时，有1 =1+hf(x,n)+...1-hf(x,n) yn+1yn+1h2即得y(x)-y =-y(x)+O(h3)n+1 n+1 2n3改进的欧拉方法用y】与y1分别表示用欧拉公式与后退的欧拉公式求得的数值解，由式(9.5)与(9.7)可知，欧拉方法与后退的欧拉方法的局部截断误差的主部只是符号相反，于是显然有， 、1， 、，一y(x)-(y+y)=O(h3)n+12n+1 n+1由此可以构造如下的计算公式yn+1F+印(X，叩+f(Xn,1，"(n=0,1,2,...N-1) (9.11)(9.11)式称为梯形方法，其局部截断误差为O(h3)，因此它是一个二阶方法，梯形方法也是隐式的，在实际应用中时常与显式的欧拉公式联合使用，构成如下的计算格式：y^=yn+hf(七,七)y。+1)=y+h[f(x,y)+f(x,y(k))]、n+1 n2nn n+1 n+1(k=0,1,2,...) (9.12)即先用欧拉方法算出初始近似值y(0)I，然后用(9.12)的第二式进行迭代，反复改进这个近似值，苴到|y(k;i)-y(%*(S为所允许的误差)为止。而把ynk+1)取作y(xn+i)的近似值yn+1，类似地计算yn+2，yn+3，…显然，如果上述迭代序列y(0)，y(1)，…收敛，其极限便满足n+1 n+1方程y“+1=yn+2[f(xn,yn)+f(xn+1，yn+1)]即序列的极限是梯形方法所得到的解y.容易证明，只要h取n+1得充分小，上述迭代过程必定收敛，事实上，将(9.11)式与(9.12)第二式相减，得y -y(k+1)=h[f(x,y)+f(x,y(k))n+1 n+1 2 n+1n+1 n+1 n+1由于f满足李普希兹条件，故有y 一y(k+1)＜一y 一y(k)n+1 n+1 2In+1 n+1因此只要h充分小，使竽＜1，迭代过程就是收敛的，但计算时需要迭代多少次，一般无法估计，故在实用时，在h取得较小的条件下，常常让计算迭代一次就结束，将其一次迭代值取作y1n+1。这时计算公式为y^=yn+hf(x,yn)v hy=y+Jf(x,y)+f(x,，(聆)]、n+1 n2nn n+1 n+1(n=0,1,2,...) (9.13)或直接写成(9.14)y=y+h［f(X,y)+f(X,y+h(x,y))］(9.14)<n+1 n2nn n+1nnny0=a(9.13)式或(9.14)式被称为改进的欧拉公式，将(9.14)式展开后与初值问题的解的泰勒展开式比较，可知其局部截断误差仍是O(h3)，即改进的欧拉方法是二阶方法，且是显式方法。通常把(9.13)式中第一式得出的初始近似值y(0)称为预测值，第二式是对n+1预测值进行一次校正，因此称这样构造的方法为预测-校正方法。［例1］用欧拉方法与改进的欧拉方法求初值问题［dy=2x<dX3另,y(0)=1在区间［0,1］上取h=0.1的数值解。［解］ 欧拉方法的计算公式为2xy+1=y+h(^英) y0=1,h=0.1n改进的欧拉方法其计算公式为2x、TOC\o"1-5"\h\zy(0)=y+h( n)n+1 n3y2nhr2x 2x.y=y+[n+ n+^]n+1 n23y2 3(y(0))2n n+1y=1,h=0.10本题的精确解为y(x)=31+X2，可用来检验数值解的精确度，列出计算结果。(请完成)。§9.2龙格库塔方法龙格-库塔(Runge-Kutta)方法，是间接利用泰勒展开的思想构造的一类数值方法。在讨论龙格一库塔方法之前，我们先介绍泰勒方法。9.2-1泰勒方法假定初值问题(9.1)，(9.2)的解山)及函数f(x,y)是足够光滑的，利用P阶泰勒公式得y(七,1)=y(「)+hy'(\)+*矿(七)+..・+^~y(p)(x)+O(hp+1)p!n=y(x)+hf(x,y(x))+务f(,y(x))+...+%f(pT)(x,y(x))+O(hp+1)当h充分小时，略去余项O饥+1)，则可导出近似值y+i的计算公式:' / \h2/ \y =y+hf(x,y)+hfix,y)n+1 nnn2!nn(9.18)' +..・+^Lf(p—1)(x,y)(9.18)p!nn算式中用到的各阶导数*)=fJ)(x,y)可以由微分方程的右端函数计算出来。例如 ^ ""

”yr=/(x,y)危«,;n)=fxG“,yn)+f(xn,yn)fy{xn,yn)y^=f,f(Xn,yn)=[f+2ff+f2f+ff\、Ln yy yJG,y)得二阶泰勒方法（9.18）式称为,阶泰勒方法。特别地，当p=1时（9.18）式就是欧拉公式（9.3）；当p=2时得二阶泰勒方法y=y=y

n+1剧L） （9.19）［例2］用泰勒方法解Iy(0)=1分别用二阶、四阶泰勒方法计算点xn=0.1,0.2,…,1.0处的数值解，并与精确解进行比较（请完成）9.2-2龙格-库塔方法的基本思想与二阶公式的推导龙格-库塔方法的基本思想是利用f（x,y）在某些点处的值的线性组合，来构造一类计算公式，使其按泰勒展开后与初值问题的解的泰勒展式比较，在尽可能多的项完全相同以确定其中的参数，从而保证算式有较高的精确度。由于避免了在算式中直接用到f（x,y）的导数，所以说龙格-库塔方法是基于间接利用泰勒展开的思想。一般龙格-库塔方法的形式为

y=y+3k+3k+...+®kn+1n11 22 rrk1=hf(xn,yn)k=hf(x+ah,y+Pk)(9.20)2 n2n211(9.20)TOC\o"1-5"\h\zk=hf(x +a h,y +P k+P k)3 n3n311 322...k=hf(x+ah,y+Pk+...+P k)r nrnr11 r,r-1r-1其中，a.，P厂，3.为常数。选取这些常数的原则，是要求(9.20)第一式右端在3,y)处作泰勒展开，且按h的幂次整理得y=y+vh+Lvh2+Lvh3+...n+1 n1 2!2 3!3与微分方程解的泰勒展式 1 J1 ..y(Xn+1)=y(x^)+fh+万fh2+互fh3+...有尽可能多的项重合，即要求\=f'v2=f'v3=f：；・・・通常把(9.20)式叫做r级(或r段)的计算公式。如果v=f(-)，对j=1,2,...,m成立，而对j=m+1时不成立，则所得的公具称为m阶的。下面我们以二阶龙格-库塔公式为例，进行具体的推导说明。设想构造如下形式的公式：y=y+3k+3k11 1122 (921)<k=hf(x”,yj I9.21)k=hf(x+ah,y+Pk)2 n2n211要求适当选取系数""a和P，使当y=心)时，(9.21)式的局部截断误差为O(h3)。为此，将k在(x,y)展开，有k=hf+h(af+Pff)+O(h2)]将上式代入(9.21)式「整理得xyny=y+h(3+3)f+h23(af+Pff)+O(h3) (9.22)n+1 n 1 2n 2xyn在解y(x)的泰勒展式y(x)=y(x)+hy'(x)+结y,f(x)+O(h3)n+1 n n2 n中，y(x)=f(x,y(x))'y'Xx)=f(x,y)+f(x,y)•f(x,y)=f+ff，即有x y xyy(x)=y(x)+hf+兽(f+ff)+O(h3) (9.23)n+1 nn2xyn将(9.22)式与(9.23)式比较，当y=y(x)时，只需取f⑦+⑦=1<响=1 (9.24)/2［晚=2则(9.21)式的局部截断误差为o(h3)。(9.24)式中包含有四个未知数，只有三个方程，其中有一个未知数可以任意取值，且显然p=a。(9.24)式的每一组解都使(9.21)式的局部截断误差为o(加)，即都使(9.21)式成为一个二阶方法。这些方法统称为二阶龙格-库塔方法。较常用的方法有两个，一个是取以=°=1,气=气=2，这就是改进的欧拉方法(9.14)。另一个是取以=b=2,①=0,①=1，得(y=y+k,\kn=hf(x,y2) (9.25)1 nnk=hf(x+h,y+与)、2n2n2(9.25)式通常称为中心点公式，二阶龙格-库塔方法每一步需要两次计算函数f的值。高阶龙格-库塔方法的推导过程很繁，这里不再一一论述。9.2-3四阶龙格-库塔方法适当选取四个点处函数f的值作线性组合，可以构造出四阶的计算公式。这些公式也有无穷多个，下