参考: パソコンで見る天体の動き, pp.131-142
function extrapolation(funcs, t0, inits, h) local unpack = unpack or table.unpack local t0 = t0 local inits = inits local function list2oneColMat(list) local oneColMat = {} for i = 1, #list do oneColMat[i] = {} oneColMat[i][1] = list[i] end return oneColMat end local 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 local function midPoints(funcs, t0, inits, h) local S = 7 -- テキストどおり 7 通り計算する。 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 local function neville(listOfLists) local numOfLists = #listOfLists local numOfHs = #listOfLists[1] local X = {} local extrapo = {} for n = 1, numOfLists do X[n] = list2oneColMat(listOfLists[n]) -- 既知の点列を 1 列行列に変換する。 for i = 2, numOfHs do for j = 2, i do X[n][i][j] = X[n][i][j-1] + (X[n][i][j-1] - X[n][i-1][j-1]) / (4^(j-1) - 1) --順次 1/2 刻みで計算する場合の公式。 end end extrapo[n] = X[n][numOfHs][numOfHs] -- 補外値だけ返す。 end return extrapo end return function() inits = neville(midPoints(funcs, t0, inits, h)) t0 = t0 + h return t0, inits end end -- 確かめる。 do -- この聯立微分方程式を補外法で解く。 local function xDot(t, x, y) return y end local function yDot(t, x, y) return t - x end -- 1 ステップだけ計算する函数を作って、 local a = extrapolation({xDot, yDot}, 0, {0, 0}, 1) -- 初期値を更新しながら必要な回数だけ繰り返す。 for i = 1, 5 do local t0, inits = a() print(t0, table.concat(inits, ", ")) end end
真値は下のとおりである。精度が高い。