参考: パソコンで見る天体の動き, p.26-29
function rkClassic(funcs, t0, inits, h, numOfDiv) local unpack = unpack or table.unpack local dim = #funcs local numOfDiv = numOfDiv or 1 local h = h / numOfDiv local hDiv2 = h / 2 local hDiv6 = h / 6 local f1, f2, f3, f4 = {}, {}, {}, {} local temps1, temps2, temps3 = {}, {}, {} local t0 = t0 local inits = inits local halfPoint = t0 + hDiv2 local nextPoint = t0 + h return function() for i = 1, numOfDiv do for j = 1, dim do f1[j] = funcs[j](t0 , unpack(inits)) ; temps1[j] = inits[j] + hDiv2 * f1[j] end for j = 1, dim do f2[j] = funcs[j](halfPoint, unpack(temps1)) ; temps2[j] = inits[j] + hDiv2 * f2[j] end for j = 1, dim do f3[j] = funcs[j](halfPoint, unpack(temps2)) ; temps3[j] = inits[j] + h * f3[j] end for j = 1, dim do f4[j] = funcs[j](nextPoint, unpack(temps3)) end for j = 1, dim do inits[j] = inits[j] + hDiv6 * (f1[j] + 2 * (f2[j] + f3[j]) + f4[j]) end t0 = nextPoint halfPoint = t0 + hDiv2 nextPoint = t0 + h end return t0, unpack(inits) end end --- 確かめる。 do -- この聯立微分方程式を解く。 local function xDot(t, x, y) return y end -- x'(t) = y(t) local function yDot(t, x, y) return t - x end -- y'(t) = t - x(t) -- 1 ステップだけ計算する函数を作って、 local rk = rkClassic({xDot, yDot}, 0, {0, 0}, 0.25, 1) -- (函数リスト, t0, 初期値リスト, 刻み幅 [, 内部分割数]) -- その函数を必要な回数だけ繰り返す。 for i = 1, 7.5/0.25 do print(rk()) end end
Nspire で検算する。2 行目が真値、3 行目が rk23() による近似値。