extrapolation.lua

extrapolation.lua

-- Lua による常微分方程式の数値解法 / 補外法
function extrapolation(funcs, t0, inits, h, numOfDiv)
   local unpack = unpack or table.unpack
   
   local numOfDiv = numOfDiv or 1
   
   local step  = h
   local h     = h / numOfDiv
   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()
             for i = 1, numOfDiv do
                inits = neville(midPoints(funcs, t0, inits, h))
                t0    = t0 + h
             end
             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