常微分方程（組）的數值解

03-14

閱讀原文

預備知識　彈簧振子受迫運動的簡單數值計算，天體運動的簡單數值計算

　　數值解微分方程或微分方程組時，一般需要先把微分方程（組）化成一階微分方程組．例如在「彈簧振子受迫運動的簡單數值計算」中列出的微分方程為

$my = alpha y - ky + f(t)$ 　　(1)

新增變數，令 $v=y′$ ，則可變為一階微分方程組

$left{egin{aligned} &v = [-alpha u - ky + f(t)]/m\ &y = v end{aligned} ight.$ 　　(2)

方程組中， $t$ 是自變數， $y$ 和 $u$ 是因變數．給出某時刻的 $y(t),u(t),t$ 就可以通過方程組求出 $u′(t)$ 和 $y′(t)$ 的數值．

　　又例如在「天體運動的簡單數值計算」中，列出的二階微分方程組為

$left{egin{aligned} &x = -frac{GM}{(x^2 + y^2)^{3/2}} x\ &y = -frac{GM}{(x^2 + y^2)^{3/2}} y end{aligned} ight.$ 　　(3)

　　新增變數，令 $v_x = x, , v_y = y$ ，上式也可變為一階微分方程組

該式中同樣 t 是自變數，其他都是因變數，給出某時刻的 $x, y, v_x, v_y, t$ 就可以由該方程組求出因變數的一階導數．

$left{egin{aligned} &x = v_x\ &v_x = -GMx/(x^2 + y^2)^{3/2}\ &y = v_y\ &v_y = -GMy/(x^2 + y^2)^{3/2} end{aligned} ight.$ 　　(4)

　　在以上兩個詞條中，我們使用了一種較為原始的方法（微分近似）．這種方法相當於把某時刻 $t$ 的各因變數代入一階微分方程組，得到 $t$ 時刻各因變數的一階導數，再通過微分近似由這些一階導數來計算其 $[t,t+Δt]$ 時間內的增量，得到 $t+Δt$ 時刻的各因變數，再代入一階方程組得到 $t+Δt$ 時刻各因變數的一階導數，如此一直循環，得到各因變數每隔 $Δt$ 時間的值．這種方法叫做歐拉法．

（剩下部分見頂部的「閱讀原文」）