RK4 与欧拉方法比较

Question

我选修了计算物理学课程，并接受了一项求解简谐振子方程 x_dot_dot = -x 的作业。我在 Julia 中编写了这两种方法的脚本

function RK4(x0 , F, T , dt) 

F1(f) = x -> f(x...)
F2(f) = x -> f((x .+ 0.5*dt*F1(f)(x))...)
F3(f) = x -> f((x .+ 0.5*dt*F2(f)(x))...)
F4(f) = x -> f((x .+ dt*F3(f)(x))...)

t = 0:dt:T
x = zeros(length(t), length(x0))
x[1,:] = x0
for n in 1:length(t) - 1
    act(F)  = map(f -> f(x[n,:]) , F)'
    x[n+1,:]= x[n,:] .+ (dt/6)*(act(F1.(F)) .+ 2*act(F2.(F)) + 2*act(F3.(F)) + act(F4.(F)))
end
return (x,t)
end

function euler(x0 , F, T , dt) 

F1(f) = x -> f(x...)

t = 0:dt:T
x = zeros(length(t), length(x0))
x[1,:] = x0
for n in 1:length(t)-1
    act(F)  = map(f -> f(x[n,:]) , F)'
    x[n+1,:]= x[n,:] .+ dt.* act(F1.(F))
end
return (x,t)
end

代码有效，并且它解决了两种方法的 ODE，当我绘制解决方案时，RK4 在以相同的时间步长尝试这两种方法时似乎没有补偿欧拉方法。当我尝试使用 RK4 时，例如比 Euler 一半的时间步长，Euler 总是显示更准确的结果，当我尝试对两种方法使用

dt = 0.01

的时间步长时，图看起来几乎完全相同。

其余代码：

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


x0 = [1,0]
f1(x,v)= v
f2(x,v) = -x

plt = plot(xlabel = L"t")

xe , t = euler(x0 , [f1 f2] , 10*2π , 0.0005)
plot!(t, xe[:,1] , label = "Euler")
xr , t = RK4(x0 , [f1 f2] , 10*2π , 0.0005)
plot!(t, xr[:,1] , label = "RK4")

# Analytical solution
#plot!(t, cos.(t))

display(plt)

谢谢！

Answer 1

RK4

不太准确的原因是因为一个错误。问题中的实现确实很混乱。甚至是非缩进的。肯定有更好、更有效的方法来实现这一点，并且 Julia 中肯定有这样做的包。

固定的

RK4

可以是：

function RK4(x0 , F, T , dt) 
    F1(f) = x -> f(x...)
    F2(f) = x -> f((x .+ 0.5*dt*[F1(g)(x) for g in F]')...)
    F3(f) = x -> f((x .+ 0.5*dt*[F2(g)(x) for g in F]')...)
    F4(f) = x -> f((x .+ dt*[F3(g)(x) for g in F]')...)

    t = 0:dt:T
    x = zeros(length(t), length(x0))
    x[1,:] = x0
    for n in 1:length(t) - 1
        act(F)  = map(f -> f(x[n,:]) , F)'
        x[n+1,:]= x[n,:] .+ (dt/6)*(act(F1.(F)) .+ 2*act(F2.(F)) + 2*act(F3.(F)) + act(F4.(F)))
    end
    return (x,t)
end

请注意

F2

、

F3

、

F4

的定义已更改。