实施 Runge Kutta 求解阻尼摆的问题

Question

我是一名高中生，正在研究 "home project" 通过使用 Runge Kutta 方法求解微分方程来制作阻尼摆的动画。

（可以在这里看到方程式：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

我已经尝试解决这个问题一段时间了，但我已经力不从心了，所以我求助于他们。感谢任何能做到的人！

Answer 1

如果您将微分方程与 RK4 实现分开，它会有所帮助。然后你可以像文档中那样实现 RK4

k1 = f1( y1, y2, x)
l1 = f2( y1, y2, x)
k2 = f1( y1 + 0.5*h*k1, y2 + 0.5*h*l1, x + 0.5*h)
l2 = f2( y1 + 0.5*h*k1, y2 + 0.5*h*l1, x + 0.5*h)
k3 = f1( y1 + 0.5*h*k2, y2 + 0.5*h*l2, x + 0.5*h)
l3 = f2( y1 + 0.5*h*k2, y2 + 0.5*h*l2, x + 0.5*h)
k4 = f1( y1 + h*k3, y2 + h*l3, x + h)
l4 = f2( y1 + h*k3, y2 + h*l3, x + h)

y1 = y1 + h/6*(k1+2*(k2+k3)+k4)
y2 = y2 + h/6*(l1+2*(l2+l3)+l4)
x = x + h

不仅有助于避免乘以h的重复，也有助于避免基点值Ome和The相加的重复。

实施 Runge Kutta 求解阻尼摆的问题

Problems implementing Runge Kutta to solve a Damped Pendulum

vb.net

simulation

runge-kutta