Runge-Kutta 4 解错误

Error in Runge-Kutta 4 solution

我正在尝试使用 Runge-Kutta 4 方法 (RK4) 求解 ODE 系统。我正在对下面的算法进行代码测试,发现解不等于解析解(而且误差很大)。下面我包含了我对 I.V.P 的代码测试。 dy/dt = f(t,y)。我曾尝试在此代码中查找错误,但无法发现它们。任何帮助深表感谢。 全局变量 [吨 dt 增长率]

turtles-own [ state ]

to setup
clear-all
create-turtles 1 [ set state 1]
set dt .01
set growth-rate .05
reset-ticks
end

to go
tick
set t t + dt
ask turtles [ set state rk4 t state dt ]    ;integrate the diff eq.
end

;differential equation to be integrated using rk4
to-report df [ t_n state_n ]   ; i.v.p. y(dot) = f(t_n, y_n)
report growth-rate * (state_n)
end


;;;;;;;function calls

to-report rk4 [ t_n state_n h ]
let k1 df t_n state_n
let k2 df (t_n + 0.5 * h) (state_n + ((h / 2) * k1))
let k3 df (t_n + 0.5 * h) (state_n + ((h / 2) * k2))
let k4 df (t_n +       h) (state_n +            k3)
let state_n+1 state_n + ((h / 6) * (k1 + (2 * k2) + (2 * k3) + k4))
report state_n+1
end

将此函数积分到 t=100,我的误差 > 7(数值解 ~156,解析解 ~148)

我认为实现基本上没问题,但是您对需要多少个报价的解释可能不正确。如果您将 go 语句更改为以下代码,它就可以工作

to go
ask turtles [ set state rk4 t state dt ]    ;integrate the diff eq.
set t t + dt
tick
if ticks = steps / dt [ stop ]
end

你应该在状态更新后 set t t + dt,因为 state_n+1 是从 state_n 在 t_n 计算的,并且首先更新时间使其基于t_n+1。然而,在实践中,这并不能解决问题(或对值产生任何真正的影响)。但是想想从 t_0 到 t_1。您需要经过 1/dt 刻度。

所以我认为当你用积分到 t=100 的例子解释问题时,你实际上是在积分到 t=101。但我不确定,因为你没有提供你模型的那部分。