nystroem.lua

nystroem.lua

-- Lua で常微分方程式を解く / 特殊な方程式に対するナイストレム法
function nystroem(funcs, t0, inits, initsDot, h, numOfDiv)
   local unpack = unpack or table.unpack  

   local t0       = t0
   local inits    = inits
   local initsDot = initsDot

   local dim      = #funcs
   local numOfDiv = numOfDiv or 1
   local h        = h / numOfDiv

   local f   = {{}, {}, {}, {}}
   local tmp = {{}, {}, {}}

   return function()
             for i = 1, numOfDiv do   
                for j = 1, dim do f[1][j] = funcs[j](unpack(inits))  ; tmp[1][j] = inits[j] + h*2/5 * initsDot[j] + (2/25) * h * h *  f[1][j]          end
                for j = 1, dim do f[2][j] = funcs[j](unpack(tmp[1])) ; tmp[2][j] = inits[j] + h*2/3 * initsDot[j] + (2/9)  * h * h *  f[1][j]          end
                for j = 1, dim do f[3][j] = funcs[j](unpack(tmp[2])) ; tmp[3][j] = inits[j] + h*4/5 * initsDot[j] + (4/25) * h * h * (f[1][j] + f[2][j]) end
                for j = 1, dim do f[4][j] = funcs[j](unpack(tmp[3]))                                                                                end

                for j = 1, dim do    inits[j] =  inits[j] + h * (initsDot[j] + (h/192) * (23 * f[1][j] + 75  * f[2][j] - 27 * f[3][j] + 25  * f[4][j])) end
                for j = 1, dim do initsDot[j] =                  initsDot[j] + (h/192) * (23 * f[1][j] + 125 * f[2][j] - 81 * f[3][j] + 125 * f[4][j])  end

                t0 = t0 + h
             end
             return t0, inits, initsDot -- 更新した t0, {更新した初期値}, {更新した 1 次微分の初期値} という形で返す。
             
          end
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}

-- 1 ステップだけ計算する函数を作って、
local a = nystroem({xDotDot, yDotDot}, 0, inits, initsDot, 0.25, _)

-- 初期値を更新しながら必要な回数だけ繰り返す。
for i = 1, 28 do
   local t0, inits = a()
   print(t0, table.concat(inits, ", "), table.concat(initsDot, ", "))
end

end