参考: パソコンで見る天体の動き, pp.131-142
function newMat(numRows, numCols, val) local mat = {} for i = 1, numRows do mat[i] = {} for j = 1, numCols do mat[i][j] = val or nil end end return mat end function midPoints(funcs, t0, inits, h, s) local temp = {} local t = {} local x = newMat(#funcs, 1, _) -- 中点法で求めた値を入れておくための 1 列行列を用意しておく。 local X = newMat(#funcs, 1, _) -- 補外にかける値を入れておくための 1 列行列を用意しておく。 for s = 1, s do local N = 2^s local hs = h/N -- 内部分割する刻み幅の設定はテキストどおり h/2^s とする。 t[1] = t0 + hs for i = 1, #funcs do x[i][1] = inits[i] + hs * funcs[i](t0, unpack(inits)) end -- これで x1, y1, ... が得られた (このスクリプトでの表現は x[1][1], x[2][1], ...)。 for n = 0, N - 1 do for i = 1, #funcs do temp[i] = x[i][n+1] end t[n+2] = (t[n] or t0) + 2 * hs for i = 1, #funcs do x[i][n+2] = (x[i][n] or inits[i]) + 2 * hs * funcs[i](t[n+1], unpack(temp)) end end -- これで {x1, x2, ..., xN+1}, {y1, y2, ..., yN+1}, ... が得られた (このスクリプトでの表現は {x[1][1], x[1][2], ..., x[1][N+1]}, {x[2][1], x[2][2], ..., x[2][N+1]})。 for i = 1, #funcs do X[i][s] = (1/4) * (x[i][N-1] + 2 * x[i][N] + x[i][N+1]) -- 補外にかける点列を計算する。 end end return X -- { {X1, X2, ...}, {Y1, Y2, ...}, ... } という形で返す (このスクリプトでの表現は { {X[1][1], X[1][2], ...}, {X[2][1], X[2][2], ...}, ... })。 end -- 確かめる。 do -- この聯立微分方程式から、補外を適用する点列を求める。 local function xDot(t, x, y) return y end local function yDot(t, x, y) return t - x end local a = midPoints({xDot, yDot}, 0, {0, 0}, 0.5, 7) for i = 1, #a do print("{"..table.concat(a[i], ", ").."}") end end