簡諧振子受迫運動的簡單數值計算

03-17

閱讀原文

預備知識　簡諧振子受迫運動，微分近似， Matlab 編程基礎

　　以下以彈簧振子的受迫運動為例，介紹一種解 $n$ 階微分方程的簡單方法，以便了解數值解微分方程的基本思想．但是這種方法誤差較大，需要大量計算才能獲得較精確的數值解．在實際運用中已有更複雜更成熟的演算法（參考 MATLAB 常微分方程（組）數值解簡介）．

　　在簡諧振子受迫運動中，列出的二階微分方程為

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

若已知初值條件（可代入任意具體數值） $dot y(0) = v_0, , y(0) = y_0$ ，且已知驅動力 $f(t)$ ，此時可以把初值條件代入式 1求出 $t = 0$ 時的加速度．

$ddot y(0) = [-alpha dot y(0) - ky(0) + f(0)]/m$ 　　(2)

接下來的一小段微小的時間 $Δt$ 內（ $Δt$ 稱為步長，步長越小誤差越小），根據微分近似，可以算出 $t=Δt$ 時刻的狀態．

$y(Delta t) = y(0) + Delta y approx y(0) + dot y(0)Delta t$ 　　(3)

$dot y(Delta t) = dot y(0) + Delta dot y approx dot y(0)+ ddot y(0)Delta t$ 　　(4)

微分近似在這裡的物理意義是在 $Δt$ 內速度和加速度都近似為常數．把 $y(Δt)$ 和 $dot y(Delta t)$ 再次代入式 1

$ddot y(Delta t) = [-alpha dot y(Delta t) - ky(Delta t) + f(Delta t)]/m$ 　　(5)

再次使用微分近似有

$y(2Delta t) = y(Delta t) + Delta y approx y(Delta t) + dot y(Delta t)Delta t$ 　　(6)

$dot y(2Delta t) = dot y(Delta t) + Delta dot y approx dot y(Delta t)+ ddot y(Delta t)Delta t$ 　　(7)

重複以上各步驟，就可以繼續得到 $y(3Δt)$ ， $y(4Δt)$ 等的近似值．在 $y-t$ 圖中把這些散點連接起來，就得到了的函數圖．

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