参考: パソコンで見る天体の動き, p.103-109
function fehlberg(funcs, t0, inits, h) local unpack = unpack or table.unpack local dim = #funcs local f1 , f2 , f3 , f4 , f5 ,f6 = {}, {}, {}, {}, {}, {} local tmp1, tmp2, tmp3, tmp4, tmp5 = {}, {}, {}, {}, {} local t0 = t0 local inits4 = {} local inits = inits return function() for j = 1, dim do f1[j] = funcs[j](t0 , unpack(inits)) ; tmp1[j] = inits[j] + h * (1/4) * f1[j] end for j = 1, dim do f2[j] = funcs[j](t0 + h * (1/4) , unpack(tmp1)) ; tmp2[j] = inits[j] + h * (1/32) * (3 * f1[j] + 9 * f2[j]) end for j = 1, dim do f3[j] = funcs[j](t0 + h * (3/8) , unpack(tmp2)) ; tmp3[j] = inits[j] + h * (1/2197) * (1932 * f1[j] - 7200 * f2[j] + 7296 * f3[j]) end for j = 1, dim do f4[j] = funcs[j](t0 + h * (12/13), unpack(tmp3)) ; tmp4[j] = inits[j] + h * ((439/216) * f1[j] - 8 * f2[j] + (3680/513) * f3[j] - (845/4104) * f4[j]) end for j = 1, dim do f5[j] = funcs[j](t0 + h , unpack(tmp4)) ; tmp5[j] = inits[j] + h * ((-8/27) * f1[j] + 2 * f2[j] - (3544/2565) * f3[j] + (1859/4104) * f4[j] - (11/40) * f5[j]) end for j = 1, dim do f6[j] = funcs[j](t0 + h * (1/2) , unpack(tmp5)) end for j = 1, dim do inits4[j] = inits[j] + h * ((25/216) * f1[j] + (1408/2565) * f3[j] + (2197/4104) * f4[j] - (1/5) * f5[j]) end -- 4 次の解 for j = 1, dim do inits [j] = inits[j] + h * ((16/135) * f1[j] + (6656/12825) * f3[j] + (28561/56430) * f4[j] - (9/50) * f5[j] + (2/55) * f6[j]) end -- 5 次の解 t0 = t0 + h return t0, inits4, inits end end -- 確かめる。 do -- この聯立微分方程式を解く。 local function xDot(t, x, y) return y end local function yDot(t, x, y) return t - x end -- この初期値で解く。 local t0 = 0 local x0, y0 = 0, 0 local h = 1/(2^2) -- 1 ステップだけ計算する函数を作って、 local fe = fehlberg({xDot, yDot}, t0, {x0, y0}, h) -- 1 ステップだけ計算する。返値は左から独立変数の値、{4 次の解}、{5 次の解} local t, a4, a5 = fe() print(t, "{"..table.concat(a4, ", ").."}", "{"..table.concat(a5, ", ").."}") end