FehlbergClass.lua

FehlbergClass.lua

-- Lua で常微分方程式を解く / フェールベルク法
Fehlberg = class()
function Fehlberg:init(funcs, t0, inits, h, tol)
   self.funcs    = funcs
   self.t0       = t0
   self.inits    = inits
   self.h        = h
   self.numOfDiv = false
   self.dim      = #funcs
   self.tol      = tol or 0.001

   self.Unpack = unpack or table.unpack
   self.Ln     = math.log
   self.Floor  = math.floor
   self.Abs    = math.abs
   self.Max    = math.max  

   self.floorB = function(num)
                    if (num > 0) and (num < 1) then 
                       return 2^self.Floor(self.Ln(num)/self.Ln(2))
                    else 
                       return self.Floor(num)
                    end
                 end
   self.hNew = function(h, err, tol)
                  if err > tol then
                     return 0.9 * h * (tol * h/(h * err))^(1/4)
                  else
                     return h
                  end
               end
   self.maxOfErr = function(listA, listB)
                      local sute = {}
                      for i = 1, #listA do
                         sute[i] = self.Abs(listA[i] - listB[i])
                      end
                      return self.Max(self.Unpack(sute))
                   end
   self.preCalc = function(numOfDiv)
                     local f   = {{}, {}, {}, {}, {}, {}}
                     local tmp = {{}, {}, {}, {}, {}}

                     local inits4   = {}
                     local preInits = {self.Unpack(self.inits)}
                     local preT0    = self.t0
                     local preH     = self.h/numOfDiv

                     for i = 1, numOfDiv do
                        for j = 1, self.dim do f[1][j] = self.funcs[j](preT0                 , self.Unpack(preInits)); tmp[1][j] = preInits[j] + preH * (1/4)    * f[1][j] end 
                        for j = 1, self.dim do f[2][j] = self.funcs[j](preT0 + preH * (1/4)  , self.Unpack(tmp[1]))  ; tmp[2][j] = preInits[j] + preH * (1/32)   * (3         * f[1][j] + 9    * f[2][j]) end
                        for j = 1, self.dim do f[3][j] = self.funcs[j](preT0 + preH * (3/8)  , self.Unpack(tmp[2]))  ; tmp[3][j] = preInits[j] + preH * (1/2197) * (1932      * f[1][j] - 7200 * f[2][j] + 7296        * f[3][j]) end
                        for j = 1, self.dim do f[4][j] = self.funcs[j](preT0 + preH * (12/13), self.Unpack(tmp[3]))  ; tmp[4][j] = preInits[j] + preH            * ((439/216) * f[1][j] - 8    * f[2][j] + (3680/513)  * f[3][j] - (845/4104)  * f[4][j]) end
                        for j = 1, self.dim do f[5][j] = self.funcs[j](preT0 + preH          , self.Unpack(tmp[4]))  ; tmp[5][j] = preInits[j] + preH            * ((-8/27)   * f[1][j] + 2    * f[2][j] - (3544/2565) * f[3][j] + (1859/4104) * f[4][j] - (11/40) * f[5][j]) end
                        for j = 1, self.dim do f[6][j] = self.funcs[j](preT0 + preH * (1/2)  , self.Unpack(tmp[5])) end

                        if numOfDiv == 1 then -- 内部分割数が 1 のときだけ 4 次の解を計算する。
                           for j = 1, self.dim do inits4[j] = preInits[j] + preH * ((25/216) * f[1][j] + (1408/2565) * f[3][j] + (2197/4104) * f[4][j] - (1/5) * f[5][j]) end
                        end
                        for j = 1, self.dim do preInits[j] = preInits[j] + preH * ((16/135) * f[1][j] + (6656/12825) * f[3][j] + (28561/56430) * f[4][j] - (9/50) * f[5][j] + (2/55) * f[6][j]) end

                        preT0 = preT0 + preH
                     end
      
                     return preT0, inits4, preInits
                  end
end
function Fehlberg:updateFE()
   self.numOfDiv   = 1
   local t, a4, a5 = self.preCalc(self.numOfDiv)                                 -- とりあえず元の刻み幅 (内部分割数 1) で 1 回だけ計算し、
   local err       = self.maxOfErr(a4, a5)  
   if err < self.tol then                                                        -- 誤差が許容範囲内であったら初期値を更新するが、
      self.t0, self.inits = t, a5
   else                                                                          -- 誤差が許容範囲外であったら、
      self.numOfDiv       = self.h/self.floorB(self.hNew(self.h, err, self.tol)) -- 新しい内部分割数を求めて、
      t, a4, a5           = self.preCalc(self.numOfDiv)                          -- その内部分割数で計算し直して初期値を更新する。
      self.t0, self.inits = t, a5 
   end
   return self
end
function Fehlberg:print()
   print(self.t0, table.concat(self.inits, ", "), self.numOfDiv)
   return self
end


--- 確かめる。
do

-- この聯立微分方程式を解く。
local function xDot(t, x, y) return y     end -- x'(t) = y(t)     
local function yDot(t, x, y) return t - x end -- y'(t) = t - x(t)

-- インスタンスを作って、   
local fe = Fehlberg({xDot, yDot}, 0, {0, 0}, 1/(2^2), 1e-9) -- (funcs, t0, inits, h, tol)

-- 必要な回数だけ更新する。   
for i = 1, 30 do
   fe:updateFE():print()
end

end