NystroemClass.lua

NystroemClass.lua

-- Lua & TI-Nspire による常微分方程式の数値解法 / ナイストレム法
Nystroem = class()
function Nystroem:init(funcs, t0, inits, initsDot, h, numOfDiv)
   self.Unpack = unpack or table.unpack  

   self.funcs    = funcs
   self.t0       = t0
   self.inits    = inits
   self.initsDot = initsDot

   self.dim      = #funcs
   self.numOfDiv = numOfDiv or 1
   self.h        = h / self.numOfDiv

   self.f   = {{}, {}, {}, {}}
   self.tmp = {{}, {}, {}}
end
function Nystroem:updateNY()
   for i = 1, self.numOfDiv do   
      for j = 1, self.dim do self.f[1][j] = self.funcs[j](self.Unpack(self.inits)) ; self.tmp[1][j] = self.inits[j] + self.h*2/5 * self.initsDot[j] + (2/25) * self.h * self.h *  self.f[1][j]                 end
      for j = 1, self.dim do self.f[2][j] = self.funcs[j](self.Unpack(self.tmp[1])); self.tmp[2][j] = self.inits[j] + self.h*2/3 * self.initsDot[j] + (2/9)  * self.h * self.h *  self.f[1][j]                 end
      for j = 1, self.dim do self.f[3][j] = self.funcs[j](self.Unpack(self.tmp[2])); self.tmp[3][j] = self.inits[j] + self.h*4/5 * self.initsDot[j] + (4/25) * self.h * self.h * (self.f[1][j] + self.f[2][j]) end
      for j = 1, self.dim do self.f[4][j] = self.funcs[j](self.Unpack(self.tmp[3])) end

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

      self.t0 = self.t0 + self.h
   end
   return self
end
function Nystroem:print()
   print(self.t0, table.concat(self.inits, ", "), table.concat(self.initsDot, ", "))
   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 ny = Nystroem({xDotDot, yDotDot}, 0, inits, initsDot, 0.25, _)

-- 必要な回数だけそのインスタンスを更新し、更新するたびに print して確認する。
for i = 1, 28 do
   ny:updateNY():print()
end

end