$B!{(B 4th-order Runge-Kutta(runge-kutta-4.rb)
# $B%k%s%2%/%C%?$G@QJ,(B
require "derivative"
def integral(time, dt, pv, pv0)
$B!!(Bif @rfirst == nil
$B!!!!(B@rfirst = false
$B!!!!(B@w = pv.dup
$B!!(Bend
$B!!(Bd1 = derivative(pv, pv0)
$B!!([email protected] = pv.spect + d1 * dt / 2.0
$B!!(Bd2 = derivative(@w, pv0)
$B!!([email protected] = pv.spect + d2 * dt / 2.0
$B!!(Bd3 = derivative(@w, pv0)
$B!!([email protected] = pv.spect + d3 * dt
$B!!(Bd4 = derivative(@w, pv0)
$B!!(Bpv.spect = pv.spect + (d1 + d2 * 2.0 + d3 * 2.0 + d4) / 6.0 * dt
$B!!(Btime = time + dt
$B!!(Breturn time, pv
end