RK4与欧拉方法求解简谐振荡器ODE的对比异常问题问询

你的RK4实现存在核心逻辑错误：将状态向量的各分量导数拆分为独立函数处理，导致中间步的状态更新完全错误。

简谐振荡器的状态向量是 $\mathbf{x} = [x, v]$，其导数向量为 $\dot{\mathbf{x}} = [v, -x]$。RK4的中间步需要用完整的导数向量来更新整个状态，而你的代码中，每个分量导数函数（f1、f2）单独计算中间状态时，错误地用单个分量的导数去更新整个状态向量。比如计算F2(f1)时，你用x .+ 0.5*dt*F1(f1)(x)，这里F1(f1)(x)是v（标量），所以状态向量会被更新为[x+0.5dt*v, v+0.5dt*v]——这完全不符合物理规律，正确的中间状态应该是[x+0.5dt*v, v+0.5dt*(-x)]。

这种错误导致你的RK4方法实际上没有正确实现四阶龙格-库塔的迭代逻辑，反而退化成了一种变形的低阶方法，自然无法体现出比欧拉法更优的精度和稳定性。

将状态向量的导数合并为一个函数，而非拆分成分量函数，让RK4方法基于完整的状态向量计算导数。以下是修正后的代码：

修正后的数值方法实现 #

function RK4(x0, F, T, dt)
    t = 0:dt:T
    x = zeros(length(t), length(x0))
    x[1, :] = x0
    for n in 1:length(t)-1
        # 计算RK4的四个增量
        k1 = dt * F(x[n, :])
        k2 = dt * F(x[n, :] + k1/2)
        k3 = dt * F(x[n, :] + k2/2)
        k4 = dt * F(x[n, :] + k3)
        # 更新状态
        x[n+1, :] = x[n, :] + (k1 + 2k2 + 2k3 + k4)/6
    end
    return (x, t)
end

function euler(x0, F, T, dt)
    t = 0:dt:T
    x = zeros(length(t), length(x0))
    x[1, :] = x0
    for n in 1:length(t)-1
        x[n+1, :] = x[n, :] + dt * F(x[n, :])
    end
    return (x, t)
end

修正后的绘图代码 #

using Plots
using LaTeXStrings
include("./methods.jl")

# 定义完整的状态导数函数
F(x) = [x[2], -x[1]]
x0 = [1, 0]

plt = plot(xlabel = L"t", title = "简谐振荡器数值解对比")

# 等步长对比：dt=0.1（故意选大一点，更易看出差异）
xe, t = euler(x0, F, 10*2π, 0.1)
plot!(t, xe[:, 1], label = "Euler (dt=0.1)", linewidth=1)
xr, t = RK4(x0, F, 10*2π, 0.1)
plot!(t, xr[:, 1], label = "RK4 (dt=0.1)", linewidth=1)

# 解析解
plot!(t, cos.(t), label = "解析解", linewidth=2, linestyle=:dash)

display(plt)

修正后：

• 等步长下，RK4的精度远高于欧拉法：欧拉法会出现明显的能量漂移（振幅随时间增大），而RK4的曲线几乎和解析解重合。
• 当RK4使用欧拉法一半的步长时，精度优势会更加显著。
• 即使步长dt=0.01，仔细对比也能发现欧拉法的微小漂移，而RK4几乎无误差。

内容的提问来源于stack exchange，提问作者Itay2924