graggClass.lua

graggClass.lua

-- Lua による微分方程式の数値解法 / グラッグ法
Gragg = class()
function Gragg:init(funcs, t0, inits, initsDot, h, numOfDiv)
   local unpack = unpack or table.unpack
   
   self.numOfDiv = numOfDiv or 1 
   
   self.step = h
   
   self.funcs = funcs
   
   self.h        = self.step / self.numOfDiv
   self.t0       = t0
   self.inits    = inits
   self.initsDot = initsDot
   
   local list2oneColMat =
   function(list)
      local oneColMat = {}
      for i = 1, #list do
         oneColMat[i] = {}
         oneColMat[i][1] = list[i]
      end
      return oneColMat
   end
   
   local newMat =
   function(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
   
   self.euler =
   function(funcs, t0, inits, initsDot, h)
      local S    = 7 -- テキストどおり 7 通り計算する。
      local temp = {}
      
      local x  = newMat(#funcs, 1, _) -- 一種のオイラー法で求めた値を入れておくための 1 列行列を用意しておく。
      local p  = newMat(#funcs, 1, _) -- 一種のオイラー法で求めた値を入れておくための 1 列行列を用意しておく。
      local p0 = {}                   -- 一種のオイラー法で求めた値を入れておくためのリストを用意しておく。
      
      local X = newMat(#funcs, 1, _) -- 補外にかける値を入れておくための 1 列行列を用意しておく。
      local U = newMat(#funcs, 1, _) -- 補外にかける値を入れておくための 1 列行列を用意しておく。
      
      for s = 1, S do
         local N  = 2^s
         local hs = h/N -- 内部分割する刻み幅の設定はテキストどおり h/2^s とする。
         
         for n = 0, N do
            for i = 1, #funcs do
               p0[i]     = initsDot[i] + (hs/2) * funcs[i](unpack(inits))
               x[i][n+1] = (x[i][n] or inits[i]) + hs * (p[i][n] or p0[i])
               temp[i]   = x[i][n+1]
            end
            for i = 1, #funcs do
               p[i][n+1] = (p[i][n] or p0[i]) + hs * funcs[i](unpack(temp))
            end
         end
         for i = 1, #funcs do
            X[i][s] = x[i][N]
            temp[i] = X[i][s]
         end
         for i = 1, #funcs do
            U[i][s]  = p[i][N] - (hs/2) * funcs[i](unpack(temp))
         end
      end
      return X, U
   end
   
   self.neville =
   function(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
end
function Gragg:updateGR()
   for i = 1, self.numOfDiv do
      local temp1, temp2        = self.euler(self.funcs, self.t0, self.inits, self.initsDot, self.h)
      self.inits, self.initsDot = self.neville(temp1), self.neville(temp2)  
      self.t0                   = self.t0 + self.h
   end
   return self
end
function Gragg:print()
   print(self.t0, table.concat(self.inits, "\t"), table.concat(self.initsDot, "\t"))   
   return self
end


-- 確かめる。
do

-- この聯立微分方程式を補外法で解く。
local function xDotDot(x, y) return -x/((math.sqrt(x^2 + y^2))^3) end
local function yDotDot(x, y) return -y/((math.sqrt(x^2 + y^2))^3) end
local inits    = {1, 0}
local initsDot = {0, 1}

-- インスタンス化して、
local gr = Gragg({xDotDot, yDotDot}, 0, inits, initsDot, 1, _)

-- 初期値を更新しながら必要な回数だけ繰り返す。
for i = 1, 5 do
   gr:updateGR()
     :print()
end

end

results.lua

1	0.54030230586814	0.8414709848079	-0.8414709848079	0.54030230586814
2	-0.41614683654714	0.90929742682568	-0.90929742682568	-0.41614683654714
3	-0.98999249660044	0.14112000805987	-0.14112000805987	-0.98999249660045
4	-0.65364362086361	-0.75680249530793	0.75680249530793	-0.65364362086361
5	0.28366218546323	-0.95892427466314	0.95892427466314	0.28366218546323