我正在尝试为高中项目使用 RK4 为阻尼摆设置动画。

描述阻尼系统的方程如下：（来自http://www.maths.tcd.ie/~smurray/Pendulumwriteup.pdf）

$$\begin{aligned} \frac{d\theta}{dt}&=\omega,\\ \frac{d\omega}{dt}&=-\beta^2\sin\theta-k\omega+A\cos\Omega. \end{aligned}$$

在哪里$B$是一个常数，$k$是阻尼系数，$A$是驱动幅度和$\Omega$是驱动频率。

我的问题是我在每次挥杆时的幅度都在增加，改变阻尼系数不会改变这一点。

我相信我的问题要么是我的编程错误，要么是我对驱动幅度和驱动频率有什么误解。
理想情况下，我希望程序能够展示简单的谐波运动（当$k = 0$) 并且还允许用户研究阻尼。

这是我第一次尝试使用 Runge Kutta 方法，我很可能错误地实现了它。

我的代码是用 Visual Basic.net 2010 编写的，如下所示：

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

Public Function y(ByRef the As Decimal, ByRef omega As Decimal, ByRef time As Decimal)
    Return omega
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 = y(The, Ome, t)
    k1 = f(The, Ome, t)
    l2 = y(The + (0.5 * h * l1), Ome + (0.5 * h * k1), t + (0.5 * h))
    k2 = f(The + (0.5 * h * l1), Ome + (0.5 * h * k1), t + (0.5 * h))
    l3 = y(The + (0.5 * h * l2), Ome + (0.5 * h * k2), t + (0.5 * h))
    k3 = f(The + (0.5 * h * l2), Ome + (0.5 * h * k2), t + (0.5 * h))
    l4 = y(The + (h * l3), Ome + (h * k3), t + h)
    k4 = f(The + (h * l3), Ome + (h * k3), 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

这是 Dis、Vel 和 Acc 图的图片（按此顺序）
如您所见，加速度在模拟过程中分崩离析（请原谅我缺乏科学术语）为什么会发生这种情况？（我还将使用我的新 RK4 实现来更新我的代码）