Visual Basic 數值分析 - Solve an arbitrary order non stiff ODE initial value problem by using a fourth-order Runge Kutta method with a constant step size.

Solve an arbitrary order , non stiff ODE initial value problem by using a fourth-order Runge Kutta method with a constant step size. (線上執行) (原始碼下載)

工程或科學上有許多問題均可以用常微分方程式 (ordinary differential equation，ODE)來表示，其中只有一小部份可以得到分析解 (analytical solution)，絕大部分必須用數值分析的技巧以求得數值解 (numerical solution)。比較常用來解 ODE 初始值問題的方法是 Euler method 和 Runge Kutta method；Euler 雖然程式簡短、跑得快，但是結果誤差大，所以並不適合高精確度的計算；一般工程上比較偏好 fourth-order Runge Kutta method，因其精確度高，且計算效率尚在可接受之範圍。Runge Kutta 一般可分固定步長(constant step size)和可變步長兩種方法，以下介紹的是使用固定步長法來解決任意階數非聯立的 ODE 初始值問題，關於 Runge Kutta method 以及 ODE 的數值分析原理請參閱市面上的書籍。

副程式：

Public Sub RungeKutta(Problem As Long, ix As Single, tx As Single, dx As Single, _
ic() As Single, ReturnValue() As Single)
     Dim i As Long, j As Long, EqNo As Long, no As Long, x As Single
     Dim f() As Single, iy() As Single, tempY() As Single, y() As Single
     EqNo = UBound(ic)
     no = CInt((tx - ix) / dx) + 1
     ReDim f(1 To EqNo, 1 To 4), iy(1 To EqNo, 1 To 3), y(1 To EqNo)
     ReDim ReturnValue(1 To EqNo, 1 To no), tempY(1 To EqNo)
     For i = 1 To EqNo
          y(i) = ic(i)
          ReturnValue(i, 1) = y(i)
     Next
     For i = 2 To no
          x = ix + (i - 2) * dx
          For j = 1 To EqNo
               If j = 1 Then f(1, 1) = fx(Problem, x, y) Else f(j, 1) = y(j - 1)
               iy(j, 1) = y(j) + (dx / 2) * f(j, 1)
               tempY(j) = iy(j, 1)
          Next
          For j = 1 To EqNo
               If j = 1 Then f(1, 2) = fx(Problem, x + (dx / 2), tempY) Else f(j, 2) = iy(j - 1, 1)
               iy(j, 2) = y(j) + (dx / 2) * f(j, 2)
               tempY(j) = iy(j, 2)
          Next
          For j = 1 To EqNo
               If j = 1 Then f(1, 3) = fx(Problem, x + (dx / 2), tempY) Else f(j, 3) = iy(j - 1, 2)
               iy(j, 3) = y(j) + dx * f(j, 3)
               tempY(j) = iy(j, 3)
          Next
          f(1, 4) = fx(Problem, x + dx, tempY)
          For j = 2 To EqNo: f(j, 4) = iy(j - 1, 3): Next
          For j = 1 To EqNo
               y(j) = y(j) + dx * (f(j, 1) / 6 + f(j, 2) / 3 + f(j, 3) / 3 + f(j, 4) / 6)
               ReturnValue(j, i) = y(j)
          Next
     Next
End Sub

Public Function fx(Problem As Long, x As Single, v() As Single) As Single
     Select Case Problem
          Case 1
               fx = 1 + v(1) / x + (v(1) / x) ^ 2
          Case 2
               fx = v(1) / x - 3 * v(2) / (x ^ 2) + 4 * v(3) / (x ^ 3) + 5 * Log(x) + 9
     End Select
End Function

參數說明：

RungeKutta 副程式中的參數必須正確地傳入，否則就算程式沒有因為錯誤而終止，其傳回的計算結果也不會正確。

Problem：代表題號。提供這個參數的好處是：不論要計算幾題 ODE 的初始值問題，都只要簡單地在 fx 函數中增加屬於該題的程式碼，完全不須更動 RungeKutta 副程式。
ix：initial x之意，代表自變數的初始值。
tx：terminated x之意，代表自變數的終值。
dx：delta x之意，代表步長大小(step size)。
ic()：initial condition 之意，是一個 1 維陣列，代表由高階 ODE 所分解出來的每一個一階ODE，其自變數=初始值時應變數之值。陣列基底為1。
ReturnValue()：即為計算結果，是一個 2 維陣列，第 1 維代表由高階 ODE 分解成數個一階ODE 的方程式編號，第二維代表該一階 ODE 相對於每個自變數的應變數之值。陣列基底為1。

fx 函數是由 RungeKutta 呼叫，必須正確地寫下能夠代表該題 ODE 的 f(x) 運算式。

Problem：同上所述。
x：自變數在計算過程中之值。
v()：是一個 1 維陣列，代表由高階 ODE 所分解出來的每一個一階 ODE，其在計算過程中的應變數之值。陣列基底為1。

應該沒有人能夠完全看懂上述的參數說明，連我自己都覺得很 confused，不妨先閱讀以下兩個簡單的範例，再來看上述的參數說明會比較清楚些。

P.S. 請到線上執行，那裡提供更方便的方式來計算ODE初始值問題。

範例一：

試解出一階 ODE 之初始值問題。

2y' = 2 + 2(y/x) + 2(y/x)², 1≦x≦3 , y(1)=0 , step size=0.2

幾乎所有 ODE 問題解決的第一步，都是將原來的 ODE 化簡成 y'=f(x)，也就是說 y' 的係數一定是1.0，將上式兩邊除 2 可得

y' = 1 + (y/x) + (y/x)², 1≦x≦3 , y(1)=0 , step size=0.2

從以上的題目敘述我們可知：Problem=1 (Problem可以任意指定)，ix=1，tx=3，dx=0.2，ic(1)=0。由於此題是一階 ODE，並不需要分解 ODE 的步驟 (如果 ODE 是二階以上就需要分解)，所以方程式的總數只有 1 條，因此在 fx 函數中 v(1) 即是代表應變數 y 在計算過程中之值，可以寫出

fx = 1 + v(1) / x + (v(1) / x) ^ 2

另外，ReturnValue(1,1)代表當 x=1 時 y 之值，ReturnValue(1,2)代表當 x=1.2 時 y 之值.........ReturnValue(1,11)代表當 x=3 時 y 之值；ReturnValue的第 1 維只有一個元素，因為方程式的總數只有1條，ReturnValue的第 2 維有11個元素，因為 x 的 step size=0.2，1~3之間總共有11個點，將以上敘述寫成程式碼如下：

Private Sub Command1_Click()
     Dim i As Long, Result() As Single
     Dim ix As Single, tx As Single, dx As Single
     Dim ic(1) As Single
     ix = 1
     tx = 3
     dx = 0.2
     ic(1) = 0
     Call RungeKutta(1, ix, tx, dx, ic, Result)
     For i = 1 To UBound(Result, 2)
          Debug.Print Result(1, i)
     Next
End Sub

執行後在VB的即時運算視窗可看見以下的數值解。

這題有分析解 y(x)=x tan(ln(x))，因此可得正確解(exactly solution)如下：

比較數值解與正確解，其誤差大約在 1×10^-6 ~1×10^-5 左右。

範例二：

試解出三階 ODE 之初始值問題。

x^³y''' - x^²y'' + 3xy' - 4y = 5x^³ln(x) + 9x^³, 1≦x≦2 , y(1)=0 , y'(1)=1 , y''(1)=3 , step size=0.1

先將上式化簡成 y'''=f(x)

y''' = y''/x - 3y'/x^² + 4y/x^³ + 5ln(x) + 9

任意的 n 階 ODE 皆可分解成 n 條聯立方程式。令 y_₁=y'， y_₂=y''，因此上述的三階 ODE 可以分解成三條聯立方程式：

第一條：y_₂' = y_₂/x - 3y_₁/x^² + 4y/x^³ + 5ln(x) + 9 , y_₂(1)=3

第二條：y_₁'=y_₂ , y_₁(1)=1

第三條：y'=y_₁ , y(1)=0

上述三條方程式的寫出順序非常重要，因為這關係到寫程式時變數值的設定。方程式的寫出順序是由 y 的最高階微分，遞減至 y 的 1 次微分。從以上三條方程式我們可知：Problem=2(Problem可以任意指定)，ix=1，tx=2，dx=0.1，ic(1)=3，ic(2)=1，ic(3)=0 (ic設定的順序即為方程式的寫出順序)；在 fx 函數中，v(1)代表第一條方程式的應變數 y_₂ 在計算過程中之值，v(2)代表第二條方程式的應變數 y_₁ 在計算過程中之值，v(3)代表第三條方程式的應變數 y 在計算過程中之值。因此可以寫出

fx = v(1) / x - 3 * v(2) / (x ^ 2) + 4 * v(3) / (x ^ 3) + 5 * Log(x) + 9

另外
ReturnValue(1,1~11)代表當x=1~2時所對應的第一條方程式的應變數 y_₂ 之值。
ReturnValue(2,1~11)代表當x=1~2時所對應的第二條方程式的應變數 y_₁ 之值。
ReturnValue(3,1~11)代表當x=1~2時所對應的第三條方程式的應變數 y 之值。

ReturnValue的第 1 維有 3 個元素，因為方程式的總數有 3 條；ReturnValue的第 2 維有 11 個元素，因為 x 的 step size=0.1，1~2之間總共有11個點，將以上敘述寫成程式碼如下：

Private Sub Command2_Click()
     Dim i As Long, Result() As Single
     Dim ix As Single, tx As Single, dx As Single
     Dim ic(1 To 3) As Single
     ix = 1
     tx = 2
     dx = 0.1
     ic(1) = 3
     ic(2) = 1
     ic(3) = 0
     Call RungeKutta(2, ix, tx, dx, ic, Result)
     For i = 1 To UBound(Result, 2)
          Debug.Print Result(3, i)
     Next
End Sub

執行後在VB的即時運算視窗可看見以下的數值解。

這題有分析解 y(x) = - x^² + x cos(ln(x)) + x sin(ln(x)) + x^³ln(x)，因此可得正確解(exactly solution)如下：

比較數值解與正確解，其誤差大約在 1×10^-6 ~1×10^-5 左右。

RungeKutta 副程式可以解決任意階數非聯立的 ODE 初始值問題，只要按照範例一以及範例二的說明，將所需條件正確地設定即可。但是這個副程式至少還有以下缺失，不過對於解決一般的問題仍可適用。

無法處理剛性(stiffness)常微分方程式。
無法處理 f(x) 有奇異點(singular point)的情形。

如果發生很嚴重的捨入誤差(roundoff error)，就必須將原來宣告成 Single 的變數改成 Double。

[ 上一個 | 首頁 | 數值分析 | 下一個 ]

Total pageview since 4/6/1999.