常微分方程数值解法在动力天文中的应用

1、 青岛科技大学 学 士 学 位 论 文 （2010年2014年）题 目： 常微分方程数值解法在动力天文中的应用学 院： 数理学院 专业班级： 信计 101 学生姓名： 潘力杨 学号：1006010129指导教师： 王斌 起讫日期： 2014年3月2014年6月 19常微分方程数值解法在动力天文中的应用前言常微分方程，作为数学分支的新发展在人类认识自然，改造自然方面发挥了巨大的力量。牛顿在研究天体力学和机械动力学时，一微分方程为工具，从而在理论上得出了行星运动的规律。后来，法国天文学家勒维烈和英国天文学家亚当斯都各自分别利用了常微分方程，演算出了当时为止的海王星位置及其运动规律。如今，常微分方程

2、更是在许多学科领域内发挥着重要作用。小行星、人造天体的运动及其规律，导弹轨道的计算，飞机空气动力学的研究，化学反应过程问题等，这些问题都可以转化为求常微分方程的解，或研究解的性质的问题。虽然常微分方程存在众多的研究和应用领域，但事实上，只有极少数诸如线性系数常微分方程和个别典型微分方程的解能用初等函数，特殊函数或它们的级数与积分表达。而在变系数常微分方程的解析求解上已经比较困难了，更别说一般的非线性常微分方程了。大多数情况下，常微分方程只能用近似方法求解，近似方法可大致分为两大类：用级数解法，逐次逼近法等求展开式的近似解析法；和在方程一些离散点上给出近似解得数值解法。在动力天文的具体应用中，事

3、实上由于所处的角度关系，我们很难宏观直接地去观测天体的具体运动状态，并直接以经典物理的力学方法去计算它的准确运动轨道及规律。由于涉及到太多的不确定因素以及运动状态的复杂性，我们只能尽量地以我们的实际观测值为参数，建立模型，运用数学的方法去研究，拟合它的运动轨迹，以此为基础才能进一步作出定性的分析。在天体运动方程的求解问题上，基本上来说也只有二体问题能给出严格解。因此为了得出满足具体观测精度的定量解，我们不得不求助于数值解法。利用数值解法求解常微分方程初值问题，一般只需要在满足规定精度的条件下求得解在若干相应点上的近似解，或是求出解的近似表达式就行，而无需使用复杂的方法试图去拟合解的表达式，或花

4、费大量的计算去得出形式解。由于它的理论简单，可操作性强且能适应各种复杂的问题条件，因此在数学建模领域有着独一无二的重要地位。在过去技术尚未得到发展的时代，在运用数值解法求解天体运动方程的过程中，处于计算精度，方程稳定性和收敛性的需要，我们往往得做出大量繁杂重复的运算从而必须花费大量的时间，可靠性也得不到保证，这也反过来制约了数值解法的应用。但现如今，随着计算机的普及以及编程运算技术的成熟运用，关于数值解法的应用领域也出现了新的契机，设定步长，循环次数的程式化运算大大缩短了我们的运算时间同时也提高了运算的精度，例如只要通过Matlab就可以完成简单方程求解问题的建模、迭代、作图、误差分析等一系列

5、的工作，大大减轻了我们的负担同时也提高了运算的效率。使得我们可以更从容地将数值解法运用到天体运动求解问题的更深邻域中去。第一章 常微分方程数值解法介绍 常微分方程可分为线性、非线性、高阶方程与方程组等类，其中线性方程包含于非线性类中，高阶方程可化为一阶方程组，若方程组中的所有未知量看作一个向量，则方程组可写成向量形式的单个方程。因此研究一阶常微分方程的初值问题 (9-1)的数值解法具有典型性，其中方程的解为。只有保证问题(9-1)的解存在惟一的前提下，其数值解法的研究才有意义。由常微分方程的基本理论，我们有：定理9.1 如果(9-1)中的满足条件（1）在区域上连续；（2）在上关于满足Lipsc

6、hitz条件，即存在常数，使得则初值问题(9-1)在区间上存在惟一的连续解。所谓数值解法，是通过常微分方程离散化而给出解在某些节点上的近似值。对于问题(9-1)，在区间引入若干节点称为由到的步长，当（常数）时称为定步长，否则称为变步长。多数情况下，采用等步长，即。记问题(9-1)的精确解为, 记的近似值为, 记为，称为问题(9-1)的数值解。 求初值问题数值解的方法一般采用步进法，即在计算出后计算，方法分为单步法和多步法。单步法是指在计算时只利用，而多步法在计算时不仅要利用，还要利用前面已经计算出来的若干个。我们称要用到和的多步方法为步方法。不论单步法和多步法，它们都有显式方法和隐式方法之分。

7、显式单步法的计算公式可以写为：，此公式右端不含。 隐式单步法的计算公式可以写为：，此公式右端含有，从而是的方程式，需要求解方程或者采用迭代法。显式多步法的计算公式可以写为。隐式多步法的计算公式可以写为 。 显式公式与隐式公式各有特点。显式公式的优点是使用方便，计算简单，效率高，其缺点是计算精度低，稳定性差。隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。 1.1 欧拉方法（显）欧拉方法的建立可以通过以下三种方法：（一）用差商近似导数，用向前差商代替，则得

8、 用其近似值代替，所得结果作为的近似值，记为，则有 这样，问题(9-1)的近似解可以表示为 (9-2)按式此式由初值经过步迭代，可逐次算出，此方程称为差分方程。式(9-2)称为求解一阶常微分方程初值问题(9-1)的欧拉公式，也称显式欧拉公式。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。（2） 用Taylor多项式近似，把在点处Taylor展开，取一次多项式近似，则得： 设步长的值较小，略去余项，并以代替，便得： （3） 将问题(9-1)中的微分方程在区间上两边积分，可得： 用,分别代替,，若对右端积分采用取左端点的矩形公式，即： 同样可得到显式欧拉公式。1.2 欧拉方法（隐）在微

9、分方程离散化时，用向后差商代替导数，即，则得到如下差分方程 (9-3)此公式称为求解问题(9-1)所得的数值解称为隐式欧拉法。隐式欧拉法与欧拉公式 (9-2)式上相似，但实际计算时却复杂得多。欧拉公式 (9-2)计算的公式中不含，但是隐式欧拉法计算的公式中含有。在求解时，为已知，是方程的根。一般说来，这是一个非线性方程，当不能精确求解得到的表达式时，我们需要运用简单迭代法来求解。迭代格式为： 由于满足Lipschitz条件，所以 由此可知，只要，迭代法就收敛到解。 1.3梯形公式 利用数值积分方法将微分方程离散化时，若用梯形公式计算式中右端积分，即 并用代替，则得计算公式 (9-4) 这就是求

10、解初值问题(9-1)的梯形公式。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 由于函数关于满足Lipschitz条件，所以 其中为Lipschitz常数。因此，当时，迭代法就收敛到解。 1.4 改进的欧拉法相比于显式欧拉法和隐式欧拉法，梯形方法提高了精度，但是梯形方法，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数的值，而迭代又要反复进行若干次，计算量很大。当函数比较复杂时，这个问题会变得更加突出。为了控制计算量，我们先用欧拉公式求得一个初步的近似值，称之为预测值，预测值的精度可能很差，再用梯形公式将它校正一次得，称为校正值。这样的预测校正系统通常称为改进的欧拉公式。即 预

11、测： 校正：为了便于编制程序上机，将上式改写成 (9-5)算法实现如表 算法9.5 输入 ，初值 计算 输出 若，置，转；否则停机。 1.5 龙格库塔（Runge-Kutta）法Runge-Kutta法由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。设初值问题(9-1)的解，由微分中值定理，我们知道，必存在，使 其中叫做在上的平均斜率。对于平均斜率，只要提供一种计算方法，该公式就给出一种数值解公式。例如，用计算，就得到一阶精度的欧拉公式；用替代，就得到隐式欧拉公式；如果用的算术平均值计算，则可得到二阶精度的梯形公式。由此可以设想，如果在上能多预测几个

12、点的斜率值，用它们的加权平均值来计算，就有望得到具有较高精度的数值公式，这就是Runge-Kutta法的基本思想。Runge-Kutta公式的一般形式是 (9-6) 其中为步长，称为方法的级数，是在点的斜率预测值，均为常数。这些常数的选取原则是使公式(9-13)具有尽可能高的精度。公式(9-13)叫做级Runge-Kutta公式，简称RK方法。 Runge-Kutta公式根据系数选择的不同，可以分为显式方法和隐式方法。当时，方法为显式方法，否则为隐式方法。显式Runge-Kutta方法可用下面的形式表示 (9-7)2级（）显式RK方法的公式为 （9-8） 这里为待定常数。 常见而二阶公式，俗称

13、中点公式为： （9-9） 类似地，对和的显式RK公式，通过更复杂的计算，可以导出三阶和四阶RK公式，其中常用的三阶和四阶RK公式分别为 (9-10)和 (9-11)式(9-22)称为经典的四阶RK公式，通常说四阶RK方法就是指用式(9-22)求解，具体算法实现如表 9.3。表 9.3算法9.2 输入 ，初值 计算输出 若，则置转 ；否则停机。 1.6 线性多步法之前的各种数值解法，都是单步法，即在每一步计算时，只要前面一个值已知的条件下就可以计算出。单步法的特点是，可以自成系统进行直接计算，因为初始条件只有一个已知，由可以计算，不必借助于其它方法，这种单步法是自开始的。如果考虑在计算值时，能够

14、比较充分地利用前面的已知信息，如和，那么就可望使所得到的更加精确。这就是多步法的基本思想。我们称要用到和的多步法为步方法。 多步法中最常用的是线性多步法，其一般公式是 (9-12)其中均为常数，。当时，上式称为显式，否则称为隐式。定义：设是初值问题(9-1)的精确解，线性多步法(9-12)在上的局部截断误差定义为： 若，则称方法(9-12)是阶方法。线性多步法的构造原理：利用Taylor展开导出线性多步公式(9-12)的基本方法是将线性多步公式(9-12)在处进行Taylor展开，然后与在处的Taylor展开式相比较，要求它们前面的项重合，由此确定参数。设初值问题(9-1)的解充分光滑，记，把

15、它在处展开成Taylor级数，则有： 及 用替代式以上两式中的，则得： 把这两个式子代入式(9-12)，得： 为了使式(9-12)具有阶精度，只须使上式的前项与在的Taylor展开式 的前项系数对应相等即可。对比关于的同次项系数，得到确定的方程组： （9-13）可见，只要式(9-12)的系数满足上式，方法(9-12)就具有阶精度。由此可得局部截断误差为： 式(9-13)共有个方程，个待定常数：，只要，即，就可以由式(9-13)，定出具有阶精度的公式(9-12)的系数. 常用的线性多部法公式有：一、阿达姆斯（Adams）公式 二、米尔恩（Milne）公式 三、海明（Hamming）公式 其中 阿

16、达姆斯（Adams）公式的构造原理为：取，有： (9-14)此时方程有9个未知数，5个方程，有四个自由变量。特别取可解得： 相应的线性多步法公式为： 因为，上式称为Adams显式公式，它是四阶公式，局部截断误差为： 如果令，由式(9-14)可得解 相应的线性多步法公式为: (9-15) 因为,式(9-15)称为Adams隐式公式，它是四阶公式，局部截断误差为: 由于线性多步法公式只有在知道了前面的之后，才能使用。所以开头的的值必须先用同阶的单步法（例如Runge-Kutta法）求出，然后才能用线性多步法公式。 第二章 常微分方程数值解法的程序实现 对于数值求解微分方程（组），MATLAB 有专

17、门的求解器可以调用。MATLAB中求微分方程的数值解命令的调用格式为：X,Y = solver(odefun,tspan ,)其中solver为MATLAB中求解微分方程的求解器：例如ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb等；odefun是要求解的显式常微分方程：；为积分区间；为初始条件。输出的X为区间上点的列向量，Y为对应于X上各个点的数值解所组成的向量。要获得问题在其他指定时间点上的解，则令（要求是单调的）。因为没有一种算法可以有效地解决所有的求解常微分方程问题，为此，MATLAB 提供了多种求解器 Solver，对于不同的常微分方程

18、问题，采用不同的Solver，一些常用的Solver总结如下：表 9.8求解器Solver求解方程类型Solver特点说明ode45非刚性 单步算法；4、5阶Runge-Kutta方程；累计截断误差达大部分场合的首选算法ode23非刚性单步算法；2、3阶Runge-Kutta方程；累计截断误差达使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到计算时间比 ode45 短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear's反向数值微分；精度中等若 ode45 失效时，可尝试使用ode23s刚性单步法；2阶 Rosebrock 算法；低

19、精度当精度较低时，计算时间比 ode15s 短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比 ode15s 短 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 MATLAB 的常用程序，其中：ode23 采用Runge-Kutta 2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度；ode45 则采用Runge-Kutta 4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度。 3.1欧拉法的Matlab程序实现Plot(x,y,x1,y1)向前欧拉法Function x,y=naeuler(dyfun,xs

20、pan,y0,h) X=xspan(1):h:xspan(2);y(1)=y0;For n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);Endx=x;y=y;x1=0:0.2:1:y1=(1+2*x1).0.5;Function y=tier(dyfun,x,y,h)y0=y;e=1e-a;K=1e+4;y=y+h*feval(dyfun,x,y);y1=y+2*e;k=1; While abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; If k>K Error(迭代发散)；

21、End End向后欧拉法（隐）Function x,y=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y(0); For n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h);Endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1).0.5;Plot(x,y,x1,y1) 隐式方法的程序去掉，整篇论文只给显式方法的程序 改进欧拉法 Function x,y=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y(0);

22、For n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2; End x=x;y=y; x1=0:0.2:1;y1=(1+2*x1).0.5; Plot(x,y,x1,y1) 2.2龙格库塔法的Matlab程序实现 三阶龙格库塔公式： function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1)

23、 = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h y(i-1)-h*K1+K2*2*h); y(i) = y(i-1)+h*(K1+4*K2+K3)/6; end format short; DELGKT3_kuta 函数运行时需要调用下列函数： function fv=Funval(f, varvec, varval) v

24、ar= findsym(f); If length(var)<4 if var(1)=varvec(1) fv=subs(f,varvec(1),varval(1); else fv=subs(f,varvec(2),varval(2); end else fv=subs(f,varvec,varval); End 四阶龙格库塔法的计算公式为： function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = f

25、indsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end format short; 同理也要先调用Funval函数（见上） 2.3adams算法的matlab程序实现 选择Adams算法中的一个显式方法，给出matlab程序实现即可第3章 动力天文中的经典问题的数值求解 3.2行星运动模型 二体问题（two-body problem）是指研究两个可以视为质点的天体在其相互之间的万有引力作用下的动力学问题。二体问题作为天体力学