Runge-Kutta 4解的误差

Question

我试图用4方法(RK4)来解决一个ODE系统。我正在对下面的算法进行代码测试，发现解决方案不等于解析解(而且误差很大)。下面我包含了I.V.P.dy/dt= f(t，y)的代码测试。我试图在这段代码中找到错误，但找不到它们。任何帮助都是非常感谢的。全球t增长率

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)。

Answer 1

我认为这个实现基本上是很好的，但是你对你需要多少个滴答的解释可能已经过时了。如果将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是从t_n的state_n中计算出来的，而时间更新首先是基于t_n+1的。然而，在实践中，这并不能解决问题(或者对值产生任何实际影响)。但是考虑一下从t_0到t_1，你需要经过1/dt的刻度。

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