v2第9章常微分方程数值解wb.docx

1、v2第9章常微分方程数值解wb第九章常微分方程数值解法在工程和科学技术的实际问题中，常需要求解常微分方程，许多实际问题的数学模型都是常微分方程或常微分方程的定解问题，如物体运动、电路振荡、化学反映及生物群体的变化等。例如：对于单摆问题(如图9.1)，选取摆长I的函数，来描述单摆运动。到下面二阶常微分方程1和质量m 1，我们希望得到摆角q的关于时间t根据力矩与角加速度的关系,经一定的简化后可得d2q dt2gsinq其中g为重力加速度。当单摆在摆动开始时刻 t t0时的初始摆角q t0qo和初始角速度凹dt t toq (t0) q0都确定时,单摆运动规律q(t)惟一确定,这就可以写成一个初值问

2、题:d2qdFgsinqq to qo, q (to) q。虽然常微分方程广泛存在于众多研究和应用领域，但是常微分方程中往往只有少数较简单和典型的微分方程(例如线性常系数常微分方程等)的解能用初等函数、特殊函数或它们的级数与积分表达 ，对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况 下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等；另一类是数值解法，它给出方程在一些离散点上的近似解。在具体求解常微分方程时，需具备某种定解条件，常微分方程和定解条件合在一起组成定解问题。定解条件有两种：一种

3、是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题；另 一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。本章主要介绍常微 分方程初值问题的基本数值方法、理论和算法。9.1基本概念9.1.1常微分方程初值问题常微分方程可分为线性、非线性、高阶方程与方程组等类，其中线性方程包含于非线性类中，高阶方 程可化为一阶方程组，若方程组中的所有未知量看作一个向量，则方程组可写成向量形式的单个方程。因此研究一阶常微分方程的初值问题dyf(x, y)，dxy(a)yo的数值解法具有典型性，其中方程的解为y : RR。a x b(9-1)只有保证问题(9-1)的解存

4、在惟一的前提下，其数值解法的研究才有意义。由常微分方程的基本理论, 我们有：定理9.1如果(9-1)中的f(x,y)满足条件(1) f (x, y)在区域 D (x, y) a x b， y 上连续；(2) f (x, y)在D上关于y满足Lipschitz条件，即存在常数 L，使得f (x, y) f (x, y) Ly y|则初值问题(9-1)在区间a,b上存在惟一的连续解 y y(x)。在本章的讨论中，我们总假定方程满足以上两个条件，从而方程总存在惟一的连续解。在实际问题的研究中，许多问题最后可以归结为一阶常微分方程组的初值问题 Yi fi (x, Y1 , Y2 丄，Ym)(i 1 ,

5、L , m)Yi(a) Yio若记 Y ( Y1, Y2,L , Ym)T , Yo (Y10, Y20,L , ymo)T, f (f1, f2,L , fm)T，则初值问题(9-2)可写成如下向量(9-2)形式的单个方程如果向量函数f(x,y)在区域D : a x b, y 使得对 x a,b, Y1, Y2 Rm，都有I f (x,yj那么问题(9-3)在a,b是存在惟一解y二y(x)。y f (x,y)y(a) y。Rm连续，且关于y满足Lipschitz条件，即存在L 0，f(x,y2) L|y1 y2，问题(9-3)与问题(9-1)形式上完全相同，故对初值问题 (9-1)所建立的各

6、种数值解法都可应用于求解问题(9-3)，只需将y换成向量y，f (x,y)换成f (x, y)即可。对于高价常微分方程的初值问题的处理，下：设有m阶常微分方程初值问题可以通过变量代换化为一阶常微分方程组初值问题，过程如(9-3)引入新变量y1y,y2y ,l ,ym(m)yy(a)y(my1y2y2y3M(m 1)、f (x, Y, Y ,L , Y )， a x b(1) (m 1) (m 1)Yo, Y (a) Yo ,L , Y (a) Yo，问题(9-4)就化为一阶常微分方程组初值问题%(a) YoY2(a) y0MYm 1 (a) y0m 2)Ym(a) y0m 1ym 1 YmYm

7、 f (x, Y1 , Y2,L,ym)(9-4)令y (Y1, Y2,l ,ym 1, ym)T, Yo (Yo, y01,l , y0m2)m 1)、TT),f (y2, yo,L , Ym, f (x,儿 y2,L ,ym)则该一阶微分方程组初值问题可写成 (9-3)的向量形式。对于问题(9-1)，在区间a,b9.1.2初值问题数值解法所谓数值解法，是通过常微分方程离散化而给出解在某些节点上的近似值。引入若干节点a xo xi X2 L Xn bhn Xn 1 Xn,( n 0,1丄,N 1)称为由Xn到Xni的步长，当hn h (常数)时称为定步长，否则称为变b a步长。多数情况下，采

8、用等步长，即 h , Xn a nh(n 1,2,L ,N)。记问题(9-1)的精确解为y(x),N记y(Xn)的近似值为yn，记f(Xn,yn)为fn， Yn 5 1,2, ,N)称为问题(9-1)的数值解。求初值问题数值解的方法一般采用步进法， 即在计算出y,i n后计算yn1，方法分为单步法和多步法。单步法是指在计算yn 1时只利用yn,而多步法在计算yn 1时不仅要利用y.,还要利用前面已经计算出来的若干个yn j, j 1,2丄，1 1。我们称要用到yn和yn j, j 1,2,L ,l 1的多步方法为I步方法。不论单步法和多步法，它们都有显式方法和隐式方法之分。显式单步法的计算公式

9、可以写为yn+1 =yn h (Xn,yn,h) (9-5)此公式右端不含yn 1。隐式单步法的计算公式可以写为yn+1二yn h (Xn, yn, yn 1,h) (9-6)在(9-6)中右端含有yn 1，从而(9-6)是yn 1的方程式，需要求解方程或者采用迭代法。显式多步法的计算公式可以写为yn+1二yn h (Xn, yn,yn-1,L ,y. l 1,h)隐式多步法的计算公式可以写为yn+1二yn h (Xn, yn 1, yn, yn-1,L , yn I 1, h), k I 1显式公式与隐式公式各有特点。显式公式的优点是使用方便，计算简单，效率高，其缺点是计算精度低，稳定性差。

10、隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再 经隐式公式迭代校正。9.2欧拉方法欧拉公式是常微分方程初值问题数值解法中最简单的方法，它的导出过程能较清楚地说明常微分方程 建立数值解法的基本思想，下面我们介绍欧拉方法。9.2.1欧拉方法欧拉方法的建立可以通过下面的三种方法：、用差商近似导数在问题(9-1)中，若用向前差商 也Q 代替y (Xn)，则得hy(Xn1)y(Xn)hf (Xn, yn(Xn)(ny(Xn)用其近似值yn代替，所得结果作为y(Xn 1)的近似值，记

11、为0,1, ,N 1)yn 1，则有yn 1 yn hf (Xn,yn) (n 0,1,L , N 1)这样，问题(9-1)的近似解可以表示为5AL 弋旳 Xj 09.2.2隐式欧拉方法和梯形方法图9.2在微分方程离散化时，如果用向后差商代替导数，即y (Xn 1)y(Xn1) y(Xn)则得到如下差分方程y 1 yn hf(Xn,yn) (n 0,1L ,N 1)(9-7) yo y(xo)按式(9-7)由初值y0经过N步迭代，可逐次算出 y1,y2, yN ,此方程称为差分方程。式 (9-7)称为求解一阶常微分方程初值问题(9-1)的欧拉公式，也称显式欧拉公式。需要说明的是，用不同的差商近

12、似导数，将得到不同的计算公式。二、用Taylor多项式近似把y(Xn 1)在Xn点处Taylor展开，取一次多项式近似，则得h2y(Xn 1) y(Xn) hy (Xn)可 yy(Xn) hf (Xn, y(Xn)h2 N y ( ) Xn,Xn 1设步长h的值较小，略去余项，并以 yn代替y(xj，便得yn 1 yn hf (Xn,yn)三、用数值积分法将问题(9-1)中的微分方程在区间xn,xn 1上两边积分，可得Xn 1y(Xn 1) y(Xn) X f(x,y(x)dx (n 0,1, , N 1) (9-8)Xn用yn 1 ,yn分别代替y(Xn 1), y(Xn)，若对右端积分采用

13、取左端点的矩形公式，即xn 1f (x,y(x)dx hf (Xn,yn)Xn同样可得到显式欧拉公式 (9-7)。类似地，对右端积分采用其它数值积分方法，又可得到不同的计算公式。以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式，其思想在研究常微分方程数值解法中具有重要意义。对于欧拉方法，从几何上看，微分方程 y f (x, y)在xoy平面上确定了一个向量场：点 (x, y)处的方向斜率为f (x, y)。如图2.2，问题(9-1)的解y y(x)代表一条过点(x0, y0)的曲线，称为积分曲线，且此 曲线上每点的切向都与向量场在这点的方向一致。从点 P0(x)

14、, y0)出发，以f(x0,y0)为斜率作一直线段，与直线x X1交于点只(为，丫1)，显然有y y hf(x,y0)，再从R出发，以彳(花，)为斜率作直线段 推进到x X2上一点P2(X2,y2)，其余类推，这样得到解曲线的一条近似曲线，它就是折线 PPP2L。因此欧拉方法又称为欧拉折线法。上面我们给出了求解一阶常微分方程初值问题 (9-1)的一种最简单的数值公式： 欧拉公式(9-7)。虽然它的精度比较低， 实践中很少采用，但它的导出过程能较清楚地说明构造数值解公式的 基本思想，且几何意义明确，因此它在理论上仍占有一定的地位。y i y hf (xn i, yn i) (n 0,1,L ,

15、N 1)yo y(x。)(9-9)此公式称为求解问题(9-1)所得的数值解称为隐式欧拉法。隐式欧拉法与欧拉公式(9-7)形式上相似，但实际计算时却复杂得多。欧拉公式 (9-7)计算不含yn 1，但是隐式欧拉法计算 yn 1的公式中含有 yn 1。在求解yn 1时，yn, Xn 1为已知,yn 1 yn hf (Xn 1, yn 1)的根。一般说来，这是一个非线性方程，当不能精确求解得到 yn 1我们需要运用简单迭代法来求解。迭代格式为yV1 yn hf (Xn,yn)ynk 11 yn hf (Xn 1, ynk1) (k 0,1,2,L )由于f (x, y)满足Lipschitz条件，所以

16、yn 1 yn 1 h f (Xn 1, yn 1) f (Xn 1, Yn 1) hL Yn 1 Yn 1由此可知，只要0 hL 1，迭代法就收敛到解 yn 1。yn 1的公式中yn 1是方程的表达式时，(9-8)中右端积分，即利用数值积分方法将微分方程离散化时，若用梯形公式计算式冷1 hx f (X, y(X)dX - f (Xn, y(Xn) f(Xn 1,y(Xn 1)n 2并用yn, yn 1代替y(Xn), y(Xn 1)，则得计算公式hyn 1 yn 2【f(Xn,yn) f (Xn 1,n1)(9-10)这就是求解初值问题(9-1)的梯形公式。梯形公式也是隐式格式，一般需用迭代

17、法求解，迭代公式为yn01 yn hf (Xn, yn)f (Xn,yn) f (Xn 1)(k0,1L )由于函数f (x, y)关于y满足Lipschitz条件，所以h2ynk11 yn 1hLkyn 1 yn 12yn 1。kf (Xn 1, yn 1) f(Xn 1,n 1)h I0 1时，迭代法就收敛到解29.2.3局部截断误差与方法的精度 为了刻画数值解的准确程度，引入局部截断误差与方法精度的概念。 定义9.1假设在某一步的数值解是准确的，即用某公式推算所得 yn 1，我们称其中L为Lipschitz常数。因此，当yn y(Xn)(这个假设称为局部化假设)。在此前提下,Rn 1y(

18、Xn 1) yn 1为该公式(即该方法)的局部截断误差。定义9.2如果某种方法的局部截断误差满足Rn 1 y( Xn1) yn 1 O(hp 1)则称该方法是 p阶方法，或具有 p阶精度。显然p越大，方法的精度越高。、欧拉法的截断误差假设问题的解y(x)充分光滑，且前n步计算结果是精确的，即yi y(Xi)，y(x) f(Xi,y(x)(i n)于是欧拉法的截断误差是R 1 y(Xny(Xnh2i)1)yn 1 y(Xn 1) yn hf(Xn,yn) y(Xn) hy (Xn)(Xn)O(h3)h2这里 y (xn)称为局部截断误差主项。显然2二、隐式欧拉法的截断误差计算如下：Rn 1 y(

19、Xn 1) yn 1 丫区 1)Yny(Xn 1) y(Xn) hy(Xn+Jh2 3y(Xn)+hy (Xn)+ y (Xn) 0(h )22Rn 1 0(h )，所以欧拉法是一阶方法。hf (Xn+i,yn+i)y(Xn) h(y(Xn)+hy (Xn)+2y(Xn) 0()2hy (Xn) O(h3)2所以隐式欧拉法是一阶方法，其局部截断误差主项是h27y(Xn) 三、梯形法的截断误差是Rn 1y(Xn 1) yn 1 y(Xn 1)ynh2f(Xn，yn) f(Xn1，yn1)hy(Xn h)亦)门丫区)y (Xn h)12 y (Xn) 0(h4)12所以梯形法是二阶方法，其局部截断

20、误差主项是h312y(Xn) 924改进的欧拉法通过上面的分析，我们看到，相比于显式欧拉法和隐式欧拉法，梯形方法提高了精度，但是梯形方法，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数 f (x, y)的值，而迭代又要反复进行若干次，计算量很大，当函数 f (x, y)比较复杂时，这个问题会变得更加突出。为了控制计算量，我们先用欧拉公式求得一个初步的近似值 yn 1，称之为预测值，预测值yn 1的精度可能很差，再用梯形公式将它校正一次得yn 1，称为校正值。这样的预测校正系统通常称为 改进的欧拉公式。即预测：yn 1 yn hf (Xn,yn)h 校正：yn 1 yn - f (X

21、n,yn) f (Xn 1,、n 1)为了便于编制程序上机，将上式改写成yp yn hf(Xnn)(9-11)yq yn hf (Xn h, yp)yn 1 2(yp yq)算法实现如表 9.1 o算法9.1输入 a,b, f (x, y), h，初值 y0N ba 门,n 0,x a,y y h计算Ypy hf (x,y)Yqy hf(x h, Yp)表9.11-(yp yq) y,x h x输出(X, y)若n N 1，置n 1 n，转；否则停机。F面计算改进欧拉法的截断误差。改进欧拉法可以改写成hYa 1 Ya f(Xn,yn) f (xn 1, Ya hf(Xn,yn)2展开式，注意到

22、 Yn Y(Xn)，有h hY(Xn) - f(Xn, Y(Xn) - f (Xn, y(Xn)2 22hfx(Xn, Y(Xn) hfy(Xn, y(Xn) f (Xn,y(Xn) O(h )h2 3Y(Xn) hy (Xn) Y (Xn) O(h )2在(xn, y(xn)处作 TaylorYn3Rn 1 y(Xn 1) Ya 1 O(h )所以改进欧拉法是二阶方法。例9.1用欧拉法、隐式欧拉法、梯形法和改进欧拉法解初值问题y y x 1, 0 x 1 dxy(0) 1此问题的精确解为 Y(x) e-X+x o解取步长h 0.1,对于此问题，欧拉法为9n 10yn +1 yn10100隐式

23、欧拉法为10n 11yn+1 yn11110梯形法为19n 1yn+1 yn 21105 10改进欧拉法为yn+1 0.905 yn0.0095n 0.1计算每种方法所得与精确解之间的误差，见表 9.2。表9.2Xn欧拉法隐式欧拉法梯形法改进欧拉法(Yn y(Xn)|)(y y(Xn)|)(y y(Xn)|)(yn y(Xn)|)0.00.00.00.00.00.134.8 10-34.3 10-57.5 101.6 10-40.28.7 10-37.7 10-31.4 10-42.9 10-40.3-21.2 10-21.0 101.9 10-44.0 10-40.41.4 10-21.3

24、10-22.2 10-44.8 10-40.51.6 10-21.4 10-2-42.5 10-45.5 100.6-21.7 10-21.6 102.7 10-45.9 10-40.7-21.8 102-21.7 10 22.9 10-46.2 10-40.81.9 10-21.7 10-23.0 10-46.5 10-40.9-21.9 10-21.8 103.1 10-46.6 10-41.01.9 10-21.8 10-23.1 10-46.6 10-4例9.2设位于坐标原点的甲舰向位于 x轴上的点A(1,0)处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度 o( o是常数)沿

25、平行于y轴的直线行驶，导弹的速度是 5 ,通过欧拉法和改进欧拉法，模拟出导弹运行的轨道曲线，并与其精确轨道进行比较。解 设导弹的轨迹曲线为 y y(x)，并设经过时间t，导弹位于点P(x, y)，乙舰位于点Q(1, ot)。由于导 弹始终对准乙舰。故此时直线 PQ就是导弹的运动轨迹曲线 0P在点P处的切线，即有dy t ydx 1 x亦即ot (1 x)y y又根据题意，弧 OP的长度为AQ的5倍，即y2dx5 t由此得1 X 2(1 x)y y 5 0 J y dx等式两边关于x求导得(1 x)y -1 y25结合初值条件y(0) 0, y(0) 0，利用常微分方程的知识，可求出此问题的精确

26、解5（1 x）5 12（18 12524F面我们考虑用数值算法来求解此问题。引入y1 y,y2 y，上述二阶常微分方程可以转化成一阶常微分方程组Y1y2y2,2y25（1 x）y1（0） 0y2（0） 09.3用欧拉法和改进欧拉法在区间 0,1求解此问题，选取h 0.01,方法求得的图像以及精确解的图像如图所示：图9.3由图可知，欧拉法和改进欧拉法都能比较好地模拟导弹的运行轨道，并且改进欧拉法更加精确。人物介绍莱昂哈德 欧拉(Leonhard Euler , 1707-1783)，瑞士数学家、自然科学家。 1707年4月15日出生于瑞士的巴塞尔，1783年9月18日于俄国圣彼得堡去逝。欧拉出生

27、于牧师家庭，自幼受父亲的教育， 13岁时入读巴塞尔大学， 15岁大学毕业， 16岁获得硕士学位。欧拉是 18世纪数学界最杰出的人物之一， 他不但为数学界做出巨大贡献， 更把数学推至几乎整个物理的领域。 他是数学史上最多产的数学家， 平均每年写出八百多页的论文， 还写了大量的力学、分析学、几何学、变分法等课本，其中无穷小分析引论、微分学原理、积分学原理等都成为数学中的经典著作。欧 拉对数学的研究如此广泛，因此在许多数学的分支中也可经常见到以他的名字命名的重要常数、公式和定理。9.3 龙格库塔( Runge-Kutta )法在数值分析中，龙格库塔法( Runge-Kutta )是用于求解常微分方程

28、重要的一类隐式或显式方法，这些技术由数学家卡尔 龙格和马丁 威尔海姆 库塔于1900年左右发明。本节我们来学习龙格-库塔法。9.3.1 Runge-Kutta 法的一般形式设初值问题(9-1)的解y y(x) C1a,b，由微分中值定理，我们知道，必存在 XnXm，使y(xn 1) y(xn ) hy ( )yn hf ( ,y( ) (9-12)yn hK *其中K叫做y(x)在Xn,Xn 1上的平均斜率。对于平均斜率 K，只要提供一种计算方法，公式 (9-12)就给 出一种数值解公式。例如，用K1 f (Xn,yn)计算K*，就得到一阶精度的欧拉公式； 用K2 f(Xn1,yn1) 替代K*，就得到隐式欧拉公式；如果用 K1, K2的算术平均值计算 K*，则可得到二阶精度的梯形公式。由此可以设想，如果在Xn,xn 1上能多预测几个点的斜率值， 用它们的加权平均值来计算 K*，就有望得到具有较高精度的数值公式，这就是 Runge-Kutta 法的基本思想。Runge-Kutta 公式的一般形式是sKi f(Xn cih,yn h aijKj)s (i 1,L ,s) (9-13)yn 1 yn h biKii1其中h为步长，s称为方法的级数， Ki是y y(x)在xn qh (0 ci 1)点的斜率预测值， G,b,aj均为 常数。这些常数的选取原则是使公式 (9-13)具