1、 青島科技大學(xué) 學(xué) 士 學(xué) 位 論 文 （2010年2014年）題 目： 常微分方程數(shù)值解法在動(dòng)力天文中的應(yīng)用學(xué) 院： 數(shù)理學(xué)院 專業(yè)班級(jí)： 信計(jì) 101 學(xué)生姓名： 潘力楊 學(xué)號(hào)：1006010129指導(dǎo)教師： 王斌 起訖日期： 2014年3月2014年6月 19常微分方程數(shù)值解法在動(dòng)力天文中的應(yīng)用前言常微分方程，作為數(shù)學(xué)分支的新發(fā)展在人類認(rèn)識(shí)自然，改造自然方面發(fā)揮了巨大的力量。牛頓在研究天體力學(xué)和機(jī)械動(dòng)力學(xué)時(shí)，一微分方程為工具，從而在理論上得出了行星運(yùn)動(dòng)的規(guī)律。后來，法國天文學(xué)家勒維烈和英國天文學(xué)家亞當(dāng)斯都各自分別利用了常微分方程，演算出了當(dāng)時(shí)為止的海王星位置及其運(yùn)動(dòng)規(guī)律。如今，常微分方程

2、更是在許多學(xué)科領(lǐng)域內(nèi)發(fā)揮著重要作用。小行星、人造天體的運(yùn)動(dòng)及其規(guī)律，導(dǎo)彈軌道的計(jì)算，飛機(jī)空氣動(dòng)力學(xué)的研究，化學(xué)反應(yīng)過程問題等，這些問題都可以轉(zhuǎn)化為求常微分方程的解，或研究解的性質(zhì)的問題。雖然常微分方程存在眾多的研究和應(yīng)用領(lǐng)域，但事實(shí)上，只有極少數(shù)諸如線性系數(shù)常微分方程和個(gè)別典型微分方程的解能用初等函數(shù)，特殊函數(shù)或它們的級(jí)數(shù)與積分表達(dá)。而在變系數(shù)常微分方程的解析求解上已經(jīng)比較困難了，更別說一般的非線性常微分方程了。大多數(shù)情況下，常微分方程只能用近似方法求解，近似方法可大致分為兩大類：用級(jí)數(shù)解法，逐次逼近法等求展開式的近似解析法；和在方程一些離散點(diǎn)上給出近似解得數(shù)值解法。在動(dòng)力天文的具體應(yīng)用中，事

3、實(shí)上由于所處的角度關(guān)系，我們很難宏觀直接地去觀測天體的具體運(yùn)動(dòng)狀態(tài)，并直接以經(jīng)典物理的力學(xué)方法去計(jì)算它的準(zhǔn)確運(yùn)動(dòng)軌道及規(guī)律。由于涉及到太多的不確定因素以及運(yùn)動(dòng)狀態(tài)的復(fù)雜性，我們只能盡量地以我們的實(shí)際觀測值為參數(shù)，建立模型，運(yùn)用數(shù)學(xué)的方法去研究，擬合它的運(yùn)動(dòng)軌跡，以此為基礎(chǔ)才能進(jìn)一步作出定性的分析。在天體運(yùn)動(dòng)方程的求解問題上，基本上來說也只有二體問題能給出嚴(yán)格解。因此為了得出滿足具體觀測精度的定量解，我們不得不求助于數(shù)值解法。利用數(shù)值解法求解常微分方程初值問題，一般只需要在滿足規(guī)定精度的條件下求得解在若干相應(yīng)點(diǎn)上的近似解，或是求出解的近似表達(dá)式就行，而無需使用復(fù)雜的方法試圖去擬合解的表達(dá)式，或花

4、費(fèi)大量的計(jì)算去得出形式解。由于它的理論簡單，可操作性強(qiáng)且能適應(yīng)各種復(fù)雜的問題條件，因此在數(shù)學(xué)建模領(lǐng)域有著獨(dú)一無二的重要地位。在過去技術(shù)尚未得到發(fā)展的時(shí)代，在運(yùn)用數(shù)值解法求解天體運(yùn)動(dòng)方程的過程中，處于計(jì)算精度，方程穩(wěn)定性和收斂性的需要，我們往往得做出大量繁雜重復(fù)的運(yùn)算從而必須花費(fèi)大量的時(shí)間，可靠性也得不到保證，這也反過來制約了數(shù)值解法的應(yīng)用。但現(xiàn)如今，隨著計(jì)算機(jī)的普及以及編程運(yùn)算技術(shù)的成熟運(yùn)用，關(guān)于數(shù)值解法的應(yīng)用領(lǐng)域也出現(xiàn)了新的契機(jī)，設(shè)定步長，循環(huán)次數(shù)的程式化運(yùn)算大大縮短了我們的運(yùn)算時(shí)間同時(shí)也提高了運(yùn)算的精度，例如只要通過Matlab就可以完成簡單方程求解問題的建模、迭代、作圖、誤差分析等一系列

5、的工作，大大減輕了我們的負(fù)擔(dān)同時(shí)也提高了運(yùn)算的效率。使得我們可以更從容地將數(shù)值解法運(yùn)用到天體運(yùn)動(dòng)求解問題的更深鄰域中去。第一章 常微分方程數(shù)值解法介紹 常微分方程可分為線性、非線性、高階方程與方程組等類，其中線性方程包含于非線性類中，高階方程可化為一階方程組，若方程組中的所有未知量看作一個(gè)向量，則方程組可寫成向量形式的單個(gè)方程。因此研究一階常微分方程的初值問題 (9-1)的數(shù)值解法具有典型性，其中方程的解為。只有保證問題(9-1)的解存在惟一的前提下，其數(shù)值解法的研究才有意義。由常微分方程的基本理論，我們有：定理9.1 如果(9-1)中的滿足條件（1）在區(qū)域上連續(xù)；（2）在上關(guān)于滿足Lipsc

6、hitz條件，即存在常數(shù)，使得則初值問題(9-1)在區(qū)間上存在惟一的連續(xù)解。所謂數(shù)值解法，是通過常微分方程離散化而給出解在某些節(jié)點(diǎn)上的近似值。對(duì)于問題(9-1)，在區(qū)間引入若干節(jié)點(diǎn)稱為由到的步長，當(dāng)（常數(shù)）時(shí)稱為定步長，否則稱為變步長。多數(shù)情況下，采用等步長，即。記問題(9-1)的精確解為, 記的近似值為, 記為，稱為問題(9-1)的數(shù)值解。 求初值問題數(shù)值解的方法一般采用步進(jìn)法，即在計(jì)算出后計(jì)算，方法分為單步法和多步法。單步法是指在計(jì)算時(shí)只利用，而多步法在計(jì)算時(shí)不僅要利用，還要利用前面已經(jīng)計(jì)算出來的若干個(gè)。我們稱要用到和的多步方法為步方法。不論單步法和多步法，它們都有顯式方法和隱式方法之分。

7、顯式單步法的計(jì)算公式可以寫為：，此公式右端不含。 隱式單步法的計(jì)算公式可以寫為：，此公式右端含有，從而是的方程式，需要求解方程或者采用迭代法。顯式多步法的計(jì)算公式可以寫為。隱式多步法的計(jì)算公式可以寫為 。 顯式公式與隱式公式各有特點(diǎn)。顯式公式的優(yōu)點(diǎn)是使用方便，計(jì)算簡單，效率高，其缺點(diǎn)是計(jì)算精度低，穩(wěn)定性差。隱式公式正好與它相反，它具有計(jì)算精度高，穩(wěn)定性好等優(yōu)點(diǎn)，但求解過程很復(fù)雜，一般采用迭代法。為了結(jié)合各自的優(yōu)點(diǎn)，通常將顯式公式與隱式公式配合使用，由顯式公式提供迭代初值，再經(jīng)隱式公式迭代校正。 1.1 歐拉方法（顯）歐拉方法的建立可以通過以下三種方法：（一）用差商近似導(dǎo)數(shù)，用向前差商代替，則得

8、 用其近似值代替，所得結(jié)果作為的近似值，記為，則有 這樣，問題(9-1)的近似解可以表示為 (9-2)按式此式由初值經(jīng)過步迭代，可逐次算出，此方程稱為差分方程。式(9-2)稱為求解一階常微分方程初值問題(9-1)的歐拉公式，也稱顯式歐拉公式。需要說明的是，用不同的差商近似導(dǎo)數(shù)，將得到不同的計(jì)算公式。（2） 用Taylor多項(xiàng)式近似，把在點(diǎn)處Taylor展開，取一次多項(xiàng)式近似，則得： 設(shè)步長的值較小，略去余項(xiàng)，并以代替，便得： （3） 將問題(9-1)中的微分方程在區(qū)間上兩邊積分，可得： 用,分別代替,，若對(duì)右端積分采用取左端點(diǎn)的矩形公式，即： 同樣可得到顯式歐拉公式。1.2 歐拉方法（隱）在微

9、分方程離散化時(shí)，用向后差商代替導(dǎo)數(shù)，即，則得到如下差分方程 (9-3)此公式稱為求解問題(9-1)所得的數(shù)值解稱為隱式歐拉法。隱式歐拉法與歐拉公式 (9-2)式上相似，但實(shí)際計(jì)算時(shí)卻復(fù)雜得多。歐拉公式 (9-2)計(jì)算的公式中不含，但是隱式歐拉法計(jì)算的公式中含有。在求解時(shí)，為已知，是方程的根。一般說來，這是一個(gè)非線性方程，當(dāng)不能精確求解得到的表達(dá)式時(shí)，我們需要運(yùn)用簡單迭代法來求解。迭代格式為： 由于滿足Lipschitz條件，所以 由此可知，只要，迭代法就收斂到解。 1.3梯形公式 利用數(shù)值積分方法將微分方程離散化時(shí)，若用梯形公式計(jì)算式中右端積分，即 并用代替，則得計(jì)算公式 (9-4) 這就是求

10、解初值問題(9-1)的梯形公式。梯形公式也是隱式格式，一般需用迭代法求解，迭代公式為 由于函數(shù)關(guān)于滿足Lipschitz條件，所以 其中為Lipschitz常數(shù)。因此，當(dāng)時(shí)，迭代法就收斂到解。 1.4 改進(jìn)的歐拉法相比于顯式歐拉法和隱式歐拉法，梯形方法提高了精度，但是梯形方法，在應(yīng)用迭代公式進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大。當(dāng)函數(shù)比較復(fù)雜時(shí)，這個(gè)問題會(huì)變得更加突出。為了控制計(jì)算量，我們先用歐拉公式求得一個(gè)初步的近似值，稱之為預(yù)測值，預(yù)測值的精度可能很差，再用梯形公式將它校正一次得，稱為校正值。這樣的預(yù)測校正系統(tǒng)通常稱為改進(jìn)的歐拉公式。即 預(yù)

11、測： 校正：為了便于編制程序上機(jī)，將上式改寫成 (9-5)算法實(shí)現(xiàn)如表 算法9.5 輸入 ，初值 計(jì)算 輸出 若，置，轉(zhuǎn)；否則停機(jī)。 1.5 龍格庫塔（Runge-Kutta）法Runge-Kutta法由數(shù)學(xué)家卡爾·龍格和馬丁·威爾海姆·庫塔于1900年左右發(fā)明。設(shè)初值問題(9-1)的解，由微分中值定理，我們知道，必存在，使 其中叫做在上的平均斜率。對(duì)于平均斜率，只要提供一種計(jì)算方法，該公式就給出一種數(shù)值解公式。例如，用計(jì)算，就得到一階精度的歐拉公式；用替代，就得到隱式歐拉公式；如果用的算術(shù)平均值計(jì)算，則可得到二階精度的梯形公式。由此可以設(shè)想，如果在上能多預(yù)測幾個(gè)

12、點(diǎn)的斜率值，用它們的加權(quán)平均值來計(jì)算，就有望得到具有較高精度的數(shù)值公式，這就是Runge-Kutta法的基本思想。Runge-Kutta公式的一般形式是 (9-6) 其中為步長，稱為方法的級(jí)數(shù)，是在點(diǎn)的斜率預(yù)測值，均為常數(shù)。這些常數(shù)的選取原則是使公式(9-13)具有盡可能高的精度。公式(9-13)叫做級(jí)Runge-Kutta公式，簡稱RK方法。 Runge-Kutta公式根據(jù)系數(shù)選擇的不同，可以分為顯式方法和隱式方法。當(dāng)時(shí)，方法為顯式方法，否則為隱式方法。顯式Runge-Kutta方法可用下面的形式表示 (9-7)2級(jí)（）顯式RK方法的公式為 （9-8） 這里為待定常數(shù)。 常見而二階公式，俗稱

13、中點(diǎn)公式為： （9-9） 類似地，對(duì)和的顯式RK公式，通過更復(fù)雜的計(jì)算，可以導(dǎo)出三階和四階RK公式，其中常用的三階和四階RK公式分別為 (9-10)和 (9-11)式(9-22)稱為經(jīng)典的四階RK公式，通常說四階RK方法就是指用式(9-22)求解，具體算法實(shí)現(xiàn)如表 9.3。表 9.3算法9.2 輸入 ，初值 計(jì)算輸出 若，則置轉(zhuǎn) ；否則停機(jī)。 1.6 線性多步法之前的各種數(shù)值解法，都是單步法，即在每一步計(jì)算時(shí)，只要前面一個(gè)值已知的條件下就可以計(jì)算出。單步法的特點(diǎn)是，可以自成系統(tǒng)進(jìn)行直接計(jì)算，因?yàn)槌跏紬l件只有一個(gè)已知，由可以計(jì)算，不必借助于其它方法，這種單步法是自開始的。如果考慮在計(jì)算值時(shí)，能夠

14、比較充分地利用前面的已知信息，如和，那么就可望使所得到的更加精確。這就是多步法的基本思想。我們稱要用到和的多步法為步方法。 多步法中最常用的是線性多步法，其一般公式是 (9-12)其中均為常數(shù)，。當(dāng)時(shí)，上式稱為顯式，否則稱為隱式。定義：設(shè)是初值問題(9-1)的精確解，線性多步法(9-12)在上的局部截?cái)嗾`差定義為： 若，則稱方法(9-12)是階方法。線性多步法的構(gòu)造原理：利用Taylor展開導(dǎo)出線性多步公式(9-12)的基本方法是將線性多步公式(9-12)在處進(jìn)行Taylor展開，然后與在處的Taylor展開式相比較，要求它們前面的項(xiàng)重合，由此確定參數(shù)。設(shè)初值問題(9-1)的解充分光滑，記，把

15、它在處展開成Taylor級(jí)數(shù)，則有： 及 用替代式以上兩式中的，則得： 把這兩個(gè)式子代入式(9-12)，得： 為了使式(9-12)具有階精度，只須使上式的前項(xiàng)與在的Taylor展開式 的前項(xiàng)系數(shù)對(duì)應(yīng)相等即可。對(duì)比關(guān)于的同次項(xiàng)系數(shù)，得到確定的方程組： （9-13）可見，只要式(9-12)的系數(shù)滿足上式，方法(9-12)就具有階精度。由此可得局部截?cái)嗾`差為： 式(9-13)共有個(gè)方程，個(gè)待定常數(shù)：，只要，即，就可以由式(9-13)，定出具有階精度的公式(9-12)的系數(shù). 常用的線性多部法公式有：一、阿達(dá)姆斯（Adams）公式 二、米爾恩（Milne）公式 三、海明（Hamming）公式 其中 阿

16、達(dá)姆斯（Adams）公式的構(gòu)造原理為：取，有： (9-14)此時(shí)方程有9個(gè)未知數(shù)，5個(gè)方程，有四個(gè)自由變量。特別取可解得： 相應(yīng)的線性多步法公式為： 因?yàn)椋鲜椒Q為Adams顯式公式，它是四階公式，局部截?cái)嗾`差為： 如果令，由式(9-14)可得解 相應(yīng)的線性多步法公式為: (9-15) 因?yàn)?式(9-15)稱為Adams隱式公式，它是四階公式，局部截?cái)嗾`差為: 由于線性多步法公式只有在知道了前面的之后，才能使用。所以開頭的的值必須先用同階的單步法（例如Runge-Kutta法）求出，然后才能用線性多步法公式。 第二章 常微分方程數(shù)值解法的程序?qū)崿F(xiàn) 對(duì)于數(shù)值求解微分方程（組），MATLAB 有專

17、門的求解器可以調(diào)用。MATLAB中求微分方程的數(shù)值解命令的調(diào)用格式為：X,Y = solver(odefun,tspan ,)其中solver為MATLAB中求解微分方程的求解器：例如ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb等；odefun是要求解的顯式常微分方程：；為積分區(qū)間；為初始條件。輸出的X為區(qū)間上點(diǎn)的列向量，Y為對(duì)應(yīng)于X上各個(gè)點(diǎn)的數(shù)值解所組成的向量。要獲得問題在其他指定時(shí)間點(diǎn)上的解，則令（要求是單調(diào)的）。因?yàn)闆]有一種算法可以有效地解決所有的求解常微分方程問題，為此，MATLAB 提供了多種求解器 Solver，對(duì)于不同的常微分方程

18、問題，采用不同的Solver，一些常用的Solver總結(jié)如下：表 9.8求解器Solver求解方程類型Solver特點(diǎn)說明ode45非剛性 單步算法；4、5階Runge-Kutta方程；累計(jì)截?cái)嗾`差達(dá)大部分場合的首選算法ode23非剛性單步算法；2、3階Runge-Kutta方程；累計(jì)截?cái)嗾`差達(dá)使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到計(jì)算時(shí)間比 ode45 短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gear's反向數(shù)值微分；精度中等若 ode45 失效時(shí)，可嘗試使用ode23s剛性單步法；2階 Rosebrock 算法；低

19、精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短 要特別的是：ode23、ode45 是極其常用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問題的解的 MATLAB 的常用程序，其中：ode23 采用Runge-Kutta 2階算法，用3階公式作誤差估計(jì)來調(diào)節(jié)步長，具有低等的精度；ode45 則采用Runge-Kutta 4階算法，用5階公式作誤差估計(jì)來調(diào)節(jié)步長，具有中等的精度。 3.1歐拉法的Matlab程序?qū)崿F(xiàn)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(迭代發(fā)散)；

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) 隱式方法的程序去掉，整篇論文只給顯式方法的程序 改進(jìn)歐拉法 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程序?qū)崿F(xiàn) 三階龍格庫塔公式： 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 函數(shù)運(yùn)行時(shí)需要調(diào)用下列函數(shù)： 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 四階龍格庫塔法的計(jì)算公式為： 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; 同理也要先調(diào)用Funval函數(shù)（見上） 2.3adams算法的matlab程序?qū)崿F(xiàn) 選擇Adams算法中的一個(gè)顯式方法，給出matlab程序?qū)崿F(xiàn)即可第3章 動(dòng)力天文中的經(jīng)典問題的數(shù)值求解 3.2行星運(yùn)動(dòng)模型 二體問題（two-body problem）是指研究兩個(gè)可以視為質(zhì)點(diǎn)的天體在其相互之間的萬有引力作用下的動(dòng)力學(xué)問題。二體問題作為天體力學(xué)