fehlberg.lua

fehlberg.lua

-- Lua で常微分方程式を解く / フェールベルク法
function fehlberg(funcs, t0, inits, h, tol)
   local unpack = unpack or table.unpack

   local t0    = t0
   local inits = inits

   local dim = #funcs
   local tol = tol or 0.001

   local Ln    = math.log
   local Floor = math.floor
   local Abs   = math.abs
   local Max   = math.max  

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

            local inits4   = {}
            local preInits = {unpack(inits)}
            local preT0    = t0
            local preH     = h/numOfDiv

            for i = 1, numOfDiv do
               for j = 1, dim do f[1][j] = funcs[j](preT0                 , unpack(preInits)); tmp[1][j] = preInits[j] + preH * (1/4)    * f[1][j] end 
               for j = 1, dim do f[2][j] = funcs[j](preT0 + preH * (1/4)  , unpack(tmp[1]))  ; tmp[2][j] = preInits[j] + preH * (1/32)   * (3         * f[1][j] + 9    * f[2][j]) end
               for j = 1, dim do f[3][j] = funcs[j](preT0 + preH * (3/8)  , 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, dim do f[4][j] = funcs[j](preT0 + preH * (12/13), 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, dim do f[5][j] = funcs[j](preT0 + preH          , 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, dim do f[6][j] = funcs[j](preT0 + preH * (1/2)  , unpack(tmp[5])) end

               if numOfDiv == 1 then -- 内部分割数が 1 のときだけ 4 次の解を計算する。
                  for j = 1, 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, 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
   
   return function()
             local numOfDiv = 1
             local t, a4, a5 = preCalc(numOfDiv)       -- とりあえず元の刻み幅で 1 回だけ計算し、
             local err = maxOfErr(a4, a5)  
             if err < tol then                         -- 誤差が許容範囲内だったら初期値を更新するが、
                t0, inits = t, a5
             else                                      -- 誤差が許容範囲外であったら、
                numOfDiv = h/floorB(hNew(h, err, tol)) -- 新しい内部分割数を求めて、
                t, a4, a5  = preCalc(numOfDiv)         -- その内部分割数で計算し直して初期値を更新する。
                t0, inits  = t, a5 
             end
             return t0, inits, numOfDiv -- 更新された t0, {更新された初期値}, 内部分割数  という形で返す。
          end
   
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)

-- 1 ステップだけ計算する函数を作って、
local a = fehlberg({xDot, yDot}, 0, {0, 0}, 1/(2^2), 1e-9) -- (funcs, t0, inits, h, tol)

-- その函数を必要な回数だけ繰り返す。
for i = 1, 30 do
   local t0, inits, numOfDiv = a() 
   print(t0, table.concat(inits, ", "), numOfDiv)
end

end