使用 RK4 进行数值自由落体分析

计算科学 数值分析 计算物理学 正则 龙格库塔
2021-12-01 04:23:53

我正在尝试计算身体自由落体的真实速度和时间。我在 Fortran 中编写了一个代码，我正在尝试使用 RK4 方法对其进行改进在此处输入图像描述

x=时间 y=总自由落体

紫线使用：

DO WHILE ((K1-K2) > y)
             y = y + t*((G*m*t/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y)))) 
             time = time + t
            END DO

绿线使用（RK4方法）：

DO WHILE ((D1-D2) > y)
           k1 = t*((G*m*(0)/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y))))
           k2 = (t/2)*((G*m*(t/2)/(2*((r-(y+k1))**2))) + sqrt(2*(G*m*(y+k1)/(r*(r-(y+k1))))))
           k3 = (t/2)*((G*m*(t/2)/(2*((r-(y+k2))**2))) + sqrt(2*(G*m*(y+k2)/(r*(r-(y+k2))))))
           k4 = t*((G*m*t/(2*((r-(y+k3))**2))) + sqrt(2*(G*m*(y+k3))/(r*(r-(y+k3)))))
             y = y + (1.0/6.0)*(k1+2*k2+2*k3+k4)
             time = time + t
            END DO

所有其他代码都是相同的（例如，步长），并且在没有 RK4 的情况下我得到了更好的结果。我想我没有正确理解 RK4，我写错了，但我看不到我的错误。我应该改变什么以获得更好的 RK4 效果还是我做错了？

1个回答
  1. 您正在使用t时间步长。用起来会更清楚dt。

  2. k1,... 应该是衍生物，而不是步骤。所以你应该使用

    k1 = ((G*m*t/(2*((r-y)**2))) + sqrt(2*(G*m*y)/(r*(r-y))))

    计算一阶导数。而且对于k2,... 也不要乘以时间。

  3. 然后在计算k2(and k3) 时不要使用y+k1但y + 0.5*dt*k1在公式中。RK4 中的前两个暂定步骤只有半步长。

为了使代码看起来更好，您应该编写一个计算加速度的函数，而不是粘贴相同的公式 4 次。减少视觉混乱也有助于思考。

你的方程式也是错误的，但这是一个单独的问题。通过这些更正，您至少应该能够从您的 Euler 和 RK 版本中获得相同的结果。

其它你可能感兴趣的问题
上一篇功率谱错误地产生负值 下一篇如果我用 WENO 在空间中离散 PDE 并用隐式方法在时间上离散，我是否需要在每个时间步求解非线性代数系统？