0

我是一名高中生，正在从事一个“家庭项目”，通过使用 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

我一直在尝试解决这个问题一段时间，而且我已经处于最后阶段，所以我求助于寻求帮助。感谢任何可以的人！

4

1 回答 1

0

如果您将微分方程与 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。

于 2015-02-26T21:51:57.270 回答