vb.net - ルンゲ・クッタを実装して減衰振子を解く問題

Question

私は高校生で、ルンゲ クッタ法を使用して微分方程式を解くことによって減衰振子をアニメーション化する「ホーム プロジェクト」に取り組んでいます。

(方程式はここで見ることができます: http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf )

私のコードでは、私の RK4 の実装は正しくないと知らされました。それを理解するのに苦労しました。プログラムは VB.net 2010 で書かれています。方程式を解くための私のコードは次のとおりです。

Public Sub RK4Solve(ByRef The As Decimal, ByRef Ome As Decimal, ByRef h As Decimal)
l1 = h * Ome
k1 = h * f(The, Ome, h)
l2 = h * (0.5 * l1) + Ome
k2 = f((The + (0.5 * h * k1)), (Ome + (0.5 * h * l1)), (t + (0.5 * h)))
l3 = h * (0.5 * l2) + Ome
k3 = f((The + (0.5 * h * k2)), (Ome + (0.5 * h * l2)), (t + (0.5 * h)))
l4 = h * l3 + Ome
k4 = f((The + (h * k3)), (Ome + (h * l3)), (t + h))
'Setting next step of variables
The = The + (h / 6 * (l1 + (2 * l2) + (2 * l3) + l4))
Ome = Ome + (h / 6 * (k1 + (2 * k2) + (2 * k3) + k4))
t += h
End Sub

各ステップに h を掛けすぎていることは承知しています。

私の完全なコード:

Public Class Form1

Dim l As Decimal = 1 'Length of rod (1m)
Dim g As Decimal = 9.81 'Gravity
Dim w As Decimal = 0 ' Angular Velocity
Dim initheta As Decimal = -Math.PI / 2 'Initial Theta
Dim theta As Decimal = -Math.PI / 2 'Theta (This one changes for the simulation)
Dim t As Decimal = 0 'Current time of the simulation
Dim h As Decimal = 0.01 'Time step
Dim b As Decimal = Math.Sqrt(g / l) 'Constant used in the function for dw/dt
Dim k As Decimal = 0 'Coefficient of Damping
Dim initialx = l * Math.Sin(initheta) 'Initial Amplitude of the pendulum

Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
End Sub

'Function for dw/dt
Public Function f(ByRef the As Decimal, ByRef omega As Decimal, ByRef time As Decimal)

    Return ((-b ^ 2) * Math.Sin(the)) - (k * omega) + (initheta * Math.Cos(omega * time))
End Function

Dim k1, k2, k3, k4, l1, l2, l3, l4 As Decimal 'Initialising RK4 variables


Public Sub RK4Solve(ByRef The As Decimal, ByRef Ome As Decimal, ByRef h As Decimal)

    l1 = h * Ome
    k1 = h * f(The, Ome, h)
    l2 = h * (0.5 * l1) + Ome
    k2 = f((The + (0.5 * h * k1)), (Ome + (0.5 * h * l1)), (t + (0.5 * h)))
    l3 = h * (0.5 * l2) + Ome
    k3 = f((The + (0.5 * h * k2)), (Ome + (0.5 * h * l2)), (t + (0.5 * h)))
    l4 = h * l3 + Ome
    k4 = f((The + (h * k3)), (Ome + (h * l3)), (t + h))

    'Setting next step of variables
    The = The + (h / 6 * (l1 + (2 * l2) + (2 * l3) + l4))
    Ome = Ome + (h / 6 * (k1 + (2 * k2) + (2 * k3) + k4))
    t += h
End Sub


'Timer ticking every 0.1s
'Time step is 0.01s to increase accuracy of results for testing
Private Sub Timer1_Tick(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Timer1.Tick
    ComboBox1.Items.Add(theta) 'Adding theta to a drop down box to test data
    RK4Solve(theta, w, h)
End Sub

Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
    Timer1.Enabled = False
End Sub
End Class

私はしばらくの間、これを解決しようとしてきましたが、最後の足を踏み入れているので、助けを求めることに頼っています. できる人に感謝します！