Bullet_General.gif (3885 bytes) Solve a set of simultaneous , first 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 RungeKuttaSF(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, k 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), 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 k = 1 To 3
               For j = 1 To EqNo
                    If k = 1 Then f(j, 1) = fx(Problem, x, y, j) Else f(j, k) = _
                    fx(Problem, x + (dx / 2), tempY, j)
                    If k = 3 Then iy(j) = y(j) + dx * f(j, 3) Else iy(j) = y(j) + (dx / 2) * f(j, k)
               Next
               For j = 1 To EqNo: tempY(j) = iy(j): Next
          Next
          For j = 1 To EqNo: f(j, 4) = fx(Problem, x + dx, tempY, j): 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, EqNo As Long) As Single
     Select Case Problem
          Case 1
                    fx = 1 + v(1) / x + (v(1) / x) ^ 2
          Case 2
               Select Case EqNo
                    Case 1
                         fx = 4 * v(1) - v(2) + v(3) - (x + 2) ^ 2
                    Case 2
                         fx = -v(1) + 3 * v(2) - 2 * v(3) + 2 * x ^ 2 + x + 15
                    Case 3
                         fx = v(1) - 2 * v(2) + 3 * v(3) - 3 * x ^ 2 + x - 10
               End Select
     End Select
End Function

參數說明：

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

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

fx 函數是由 RungeKuttaSF 呼叫，必須正確地寫下能夠代表聯立 ODE 中的每一個一階 ODE 的 f(x) 運算式。

  • Problem：同上所述。
  • x：自變數在計算過程中之值。
  • v()：是一個 1 維陣列，代表聯立 ODE 中的每一個一階ODE，其在計算過程中的應變數之值。陣列基底為1。
  • EqNo：代表聯立 ODE 中的方程式編號，不同的方程式必須使用不同的 f(x) 運算式。

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

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

範例一：

雖然 RungeKuttaSF 是用來解決一階聯立的 ODE 初始值問題，但是非聯立的一階 ODE 初始值問題亦可解出，我們就用上次解過的問題當作範例。

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

2y' = 2 + 2(y/x) + 2(y/x)2 ,  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)2 ,  1≦x≦3  ,   y(1)=0  ,  step size=0.2

從以上的題目敘述我們可知：Problem=1 (Problem可以任意指定)，ix=1，tx=3，dx=0.2，ic(1)=0。由於此題是非聯立一階 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 ic(1 To 3) As Single, dx As Single, ix As Single, tx As Single
     Dim i As Long, j As Long, Result() As Single
     ix = 1
     tx = 3
     dx = 0.2
     ic(1) = 0
     Call RungeKuttaSF(1, ix, tx, dx, ic, Result)
     For i = 1 To UBound(Result, 2)
          Debug.Print Result(1, i)
     Next
End Sub

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

0 
 0.2212457 
 0.4896842 
 0.8127522 
 1.199432 
 1.661265 
 2.213469 
 2.876494 
 3.678379 
 4.658506 
 5.873839

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

0
0.2212428
0.4896817
0.8127527
1.1994386
1.6612818
2.2135018
2.8765514
3.6784753
4.6586651
5.8741000

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

　

範例二：

試解出一階聯立 ODE 的初始值問題。

第一條ODE：y1' = 4y1 - y2 + y3 - (x+2)2  ,  0≦x≦1   ,  y1(0)=3  ,  step size=0.1

第二條ODE：y2' = -y1 + 3y2 - 2y3 + 2x2 + x +15 ,  0≦x≦1   ,  y2(0)= -7 

第三條ODE：y3' = y1 - 2y2 + 3y3 - 3x2 + x - 10 ,  0≦x≦1   ,  y3(0)= -2 

有一點要注意，如果 y' 的係數不是1，一定要化簡為1。上述一階聯立 ODE 方程式的寫出順序並不重要，但是一旦寫定就不要改變，因為這關係到寫程式時變數值的設定。從以上三條方程式我們可知：Problem=2 (Problem可以任意指定)，ix=0，tx=1，dx=0.1，ic(1)=3，ic(2)= -7，ic(3)= -2 (ic設定的順序即為方程式的寫出順序)；在 fx 函數中，v(1)代表第一條方程式的應變數 y1 在計算過程中之值，v(2)代表第二條方程式的應變數 y2 在計算過程中之值，v(3)代表第三條方程式的應變數 y3 在計算過程中之值。因此可以寫出

Select Case EqNo
     Case 1
          fx = 4 * v(1) - v(2) + v(3) - (x + 2) ^ 2
     Case 2
          fx = -v(1) + 3 * v(2) - 2 * v(3) + 2 * x ^ 2 + x + 15
     Case 3
          fx = v(1) - 2 * v(2) + 3 * v(3) - 3 * x ^ 2 + x - 10
End Select

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

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

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

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

3        -7    -2         
 4.621062   -7.681889    -1.728795     
 7.161582   -8.938198    -0.907413     
 11.26138   -11.28263     0.9731913    
 18.04544   -15.66935     4.862049     
 29.5083    -23.8617      12.51682      
 49.20923   -39.1062      27.17773     
 83.53179   -67.36343     54.79842     
 143.9689   -119.5546     106.2924     
 251.2729   -215.6591     201.6307     
 442.9977   -392.187      377.3139

這題有分析解
y1(x) = e6x + 2e3x + x
y2(x) = -e6x  + e3x - 2ex - 5
y3(x) = e6x - e3x - 2ex + x2
因此可得正確解(exactly solution)如下：

 3        -7    -2         
 4.621837   -7.682602    -1.728082     
 7.164355   -8.940804    -0.9048073    
 11.26885   -11.28976     0.9803271    
 18.06341   -15.68671     4.87941      
 29.54892   -23.90129     12.55641     
 49.29753   -39.19283     27.26435     
 83.71867   -67.54766     54.98265     
 144.3568   -119.9383     106.6762     
 252.0658   -216.4459     202.4174     
 444.5999   -393.7798     378.9067

比較數值解與正確解，可發現有相當程度的誤差，因為這題有輕微的 stiff 傾向，減少 step size 將可得到更準確的結果。

RungeKuttaSF 副程式可以解決一階聯立的 ODE 初始值問題，只要按照範例一以及範例二的說明，將所需條件正確地設定即可。同時 RungeKuttaSF 也能解決高階聯立的 ODE 初始值問題，只要將每個高階 ODE 都降成一階，再利用上述法則即可解決，不過我還沒碰過高階聯立的 ODE 初始值問題，如果誰有這方面的例題，請告訴我，謝謝。RungeKuttaSF 副程式至少還有以下缺失，不過對於解決一般的問題仍可適用。

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

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

　

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

This page was written by Jaric on Sep. 18, 1998. All rights reserved.

Total pageview since 4/6/1999.