一種二階常微分方程的數值解法

戶永清,陳正文

（1.四川文理學院智能制造學院,四川 達州 635000; 2.巴中市第二中學,四川 巴中 636000）

大多數微分方程,是沒有解析解,或者很難求得解析解,在科學與工程技術中求微分方程的數值解顯得很重要,二階常微分方程的數值解一般采用歐拉法求解,精度不高,本文找到一種數值解法,能較好解決這一問題.

一階常微分方程的數值解法常用有歐拉法、龍格—庫塔法,［1］尤其以四階龍格-庫塔法（Rk4）最為著名,因其高精度,易編程性而被廣泛采用.本文給出了一種二階常微分方程的數值解法,第一步:將二階微分方程降階為一階微分方程的標準形式;第二步:利用四階龍格-庫塔法求出迭代區間中點、右端點一階導數數值解;第三步:利用迭代區間左、中、右三點一階導數值結合辛普森數值積分法迭代出原函數數值解.其中第二步使用龍格—庫塔法,公式中原函數值的表示,利用泰勒展式恰好可以保留到二階項,比起只保留一階項,提高了精度.

1 一階常微分方程的Rk4法［1］

2 二階常微分方程的數值解法

設有:y=y（x）,二階常微分方程可化為y″=f（x,y′,y）的標準形式.采用降階法:令v=y′則:v′=f（x,v,y）,其中v=v（x,y）（y=y（x））,初值條件:yk=y（x0）,vk=v（x0）,yk″=v′（x0）=f（x0,v（x0）,y（x0））.這個問題就變成了求兩個一階微分方程的數值解,由于兩個微分方程的各量彼此依賴,所以并不是簡單的解兩個獨立的一階微分方程.由于y′=v=v（x,y（x））=v（x）.假設其圖象如圖1所示.

圖1 v（x）圖像

該方法的誤差分析:Rk4 的截斷誤差為O（h5）,具有四階精度,泰勒展式截取前三項的截斷誤差為O（h3）具有二階精度,辛普森積分的截斷誤差為O（h4）,具有三階精度,而一般采用歐拉法求解二階常微分方程數值解,僅具有一階精度,因而相比歐拉法該方法具有較好的精度.

3 應 用

物理學中許多運動,都可以建立起二階常微分方程,可用上述方法求得數值解,并可以和解析解比較來檢驗該方法的可靠性.

3.1 有心力作用下質點的運動［3-4］

物體所受作用力始終通過某一點,稱為有心力,該點叫作力心,有心力作用下質點運動的求解,最好建立極坐標系,以力心為極點,如圖2所示.

圖2 極坐標系中的有心力

設質點P 的位置矢量為r→,r?、θ?分別為徑向單位矢量、橫向單位矢量,則r→=rr?（r為極徑）.

3.1.1 天體在與距離二次方成反比的萬有引力作用下的運動求解

數值解軌跡如圖3,解析解［3］如圖4,初始徑向速度分量不等于0時的數值解軌道如圖5.

圖3 天體運動數值解軌跡

圖5 天體數值解軌道圖

3.1.2 天體在與距離成反比的有心引力作用下的運動仿真

天體在與距離成反比的有心引力作用下的運動,理論上可證明其軌道是周期性不閉合的,數值解顯示的軌道與理論相符,如圖6所示.

圖6 天體運動數值解軌道

3.2 彈簧振子的運動求解

微分方程的建立,以水平彈簧振子為例,如圖7,設彈簧的勁度系數為k,振子質量為m,O 為平衡位置,振子所受阻力滿足f=－kfv,有:=bx,初值條件:xk=x0,Vk=v,VB編程,主程序略.

圖7 水平彈簧振子

數值解如圖8 所示,解析解和數值解圖像對比（藍色圖像為解析解圖像,紅色圖像為數值解圖像）如圖9所示.

圖8 水平彈簧振子數值解圖像

圖9 水平彈簧振子解析解和數值解對比圖像

結束語

通過以上物理學上的三個具體的運動,建立二階常微分方程,求其數值解并與解析解比較,發現用本文的數值解法與真實解符合的很好,具有較高的精度.