Fortran 中的龙格库塔

Question

我正在尝试在 Fortran 中实现 Runge Kutta 方法，但遇到了收敛问题。我不知道应该显示多少代码，所以我会详细描述问题，请指导我应该 add/remove to/from post让它负责。

我有一个球的位置和速度的 6 维矢量，以及相应的 diff 系统。均衡器。描述运动方程，我想从中计算球的轨迹，并比较 RK 方法不同阶数的结果。

让我们关注三阶 RK。我使用的模型实现如下：

k1 = h * f(vec_old,omega,phi)
k2 = h * f(vec_old + 0.5d0 * k1,omega,phi)
k3 = h * f(vec_old + 2d0 * k2 - k1,omega,phi)
vec = vec_old + (k1 + 4d0 * k2 + k3) / 6d0

其中 f 是构成运动方程的函数（或等价于我的差分方程系统的 RHS）。请注意 f 与时间无关，因此只有 1 个参数。 h 扮演小时间步dt的角色。

如果我们希望在有限时间内计算球的轨迹 total_time，并允许总误差为 epsilon，那么我们需要确保每一步都占比例分数错误。对于第一步，我执行了以下操作：

vec1 = solve(3,vec_old,h,omega,phi)
vec2 = solve(3,vec_old,h/2d0,omega,phi)
do while (maxval((/(abs(vec1(i) - vec2(i)),i=1,6)/)) > eps * h / (tot_time - current_time))
    h = h / 2d0
    vec1 = solve(3,vec_old,h,omega,phi)
    vec2 = solve(3,vec_old,h/2d0,omega,phi)
end do
vec = (8d0/7d0) * vec2 - (1d0/7d0) * vec1

其中solve(3,vec_old,h,omega,phi)是计算上述单个RK步骤的函数。 3 表示我们正在使用的 RK 顺序，vec_old 是位置速度向量的当前状态，h, h/2d0 都表示正在使用的时间步长，omega,phi 只是f 的一些额外参数。最后，第一步我们设置 current_time = 0d0.

重点是，如果我们使用三阶 RK，我们应该在 $O(h^3)$ 中有误差，因此比 h 中的线性下降更快。因此，我们应该期望 while 循环最终会停止足够小的 h.

我的问题是循环不收敛，甚至不收敛 - 比率

maxval(...) / eps * (...)

几乎保持不变，一直到 eps * h / (tot_time - current_time)) 由于精度有限而变为零。

为了完整起见，这是我对 f 的定义：

function f(vec_old,omega,phi) result(vec)
    real(8),intent(in) :: vec_old(6),omega,phi
    real(8) :: vec(6)
    real(8) :: v,Fv

    v = sqrt(vec_old(4)**2+vec_old(5)**2+vec_old(6)**2)
    Fv = 0.0039d0 + 0.0058d0 / (1d0 + exp((v-35d0)/5d0))
    vec(1) = vec_old(4)
    vec(2) = vec_old(5)
    vec(3) = vec_old(6)
    vec(4) = -Fv * v * vec_old(4) + 4.1d-4 * omega * (vec_old(6)*sin(phi) - vec_old(5)*cos(phi))
    vec(5) = -Fv * v * vec_old(5) + 4.1d-4 * omega * vec_old(4)*cos(phi)
    vec(6) = -Fv * v * vec_old(6) - 4.1d-4 * omega * vec_old(4)*sin(phi) - 9.8d0
end function f

有人知道为什么 while 循环不收敛吗？如果需要其他任何东西（输出、其他代码片段等），请告诉我，我会添加它。此外，如果需要修剪，我会剪掉任何被认为不必要的东西。谢谢！

Answer 1

要使用半步法计算步进误差，您需要在两种情况下计算 t+h 处的近似值，这意味着步长为 h/2 的两个步长。现在，您将 t+h 处的近似值与 t+h/2 处的近似值进行比较，这会给出大小为 f(vec(t+h/2))*h/2.

的误差

因此改为三步程序

vec1 = solve(3,vec_old,h,omega,phi)
vec2 = solve(3,vec_old,h/2d0,omega,phi)
vec2 = solve(3,vec2   ,h/2d0,omega,phi)

在两个位置，vec2-vec1 的差异应该是 h^4。

Fortran 中的龙格库塔

Runge Kutta in Fortran

fortran

fortran90

runge-kutta