rkShanks8.lua

rkShanks8.lua

function rkShanks8(funcs, t0, inits, h, numOfDiv)
   local unpack = unpack or table.unpack  

   local step  = h 
   local t0    = t0
   local inits = inits

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

   local fnc = {{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}
   local tmp = {{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}

   return function()
             for i = 1, numOfDiv do   
                for j = 1, dim do fnc[01][j] = funcs[j](t0        , unpack(inits))   ; tmp[01][j] = inits[j] + (h/9)    * (        fnc[01][j]) end
                for j = 1, dim do fnc[02][j] = funcs[j](t0 + h/9  , unpack(tmp[01])) ; tmp[02][j] = inits[j] + (h/24)   * (        fnc[01][j] + 3 * fnc[02][j]) end
                for j = 1, dim do fnc[03][j] = funcs[j](t0 + h/6  , unpack(tmp[02])) ; tmp[03][j] = inits[j] + (h/16)   * (        fnc[01][j]                  + 3  * fnc[03][j]) end
                for j = 1, dim do fnc[04][j] = funcs[j](t0 + h/4  , unpack(tmp[03])) ; tmp[04][j] = inits[j] + (h/500)  * (29    * fnc[01][j]                  + 33 * fnc[03][j] - 12    * fnc[04][j]) end
                for j = 1, dim do fnc[05][j] = funcs[j](t0 + h/10 , unpack(tmp[04])) ; tmp[05][j] = inits[j] + (h/972)  * (33    * fnc[01][j]                                    + 4     * fnc[04][j] + 125   * fnc[05][j]) end
                for j = 1, dim do fnc[06][j] = funcs[j](t0 + h/6  , unpack(tmp[05])) ; tmp[06][j] = inits[j] + (h/36)   * (-21   * fnc[01][j]                                    + 76    * fnc[04][j] + 125   * fnc[05][j] - 162   * fnc[06][j]) end
                for j = 1, dim do fnc[07][j] = funcs[j](t0 + h/2  , unpack(tmp[06])) ; tmp[07][j] = inits[j] + (h/243)  * (-30   * fnc[01][j]                                    - 32    * fnc[04][j] + 125   * fnc[05][j]                      + 99   * fnc[07][j]) end
                for j = 1, dim do fnc[08][j] = funcs[j](t0 + h*2/3, unpack(tmp[07])) ; tmp[08][j] = inits[j] + (h/324)  * (1175  * fnc[01][j]                                    - 3456  * fnc[04][j] - 6250  * fnc[05][j] + 8424  * fnc[06][j] + 242  * fnc[07][j] - 27  * fnc[08][j]) end
                for j = 1, dim do fnc[09][j] = funcs[j](t0 + h/3  , unpack(tmp[08])) ; tmp[09][j] = inits[j] + (h/324)  * (293   * fnc[01][j]                                    - 852   * fnc[04][j] - 1375  * fnc[05][j] + 1836  * fnc[06][j] - 118  * fnc[07][j] + 162 * fnc[08][j] + 324  * fnc[09][j]) end
                for j = 1, dim do fnc[10][j] = funcs[j](t0 + h*5/6, unpack(tmp[09])) ; tmp[10][j] = inits[j] + (h/1620) * (1303  * fnc[01][j]                                    - 4260  * fnc[04][j] - 6875  * fnc[05][j] + 9990  * fnc[06][j] + 1030 * fnc[07][j]                                        + 162  * fnc[10][j]) end
                for j = 1, dim do fnc[11][j] = funcs[j](t0 + h*5/6, unpack(tmp[10])) ; tmp[11][j] = inits[j] + (h/4428) * (-8595 * fnc[01][j]                                    + 30720 * fnc[04][j] + 48750 * fnc[05][j] - 66096 * fnc[06][j] + 378  * fnc[07][j] - 729 * fnc[08][j] - 1944 * fnc[09][j] - 1296 * fnc[10][j] + 3240 * fnc[11][j]) end
                for j = 1, dim do fnc[12][j] = funcs[j](t0 + h    , unpack(tmp[11])) end

                for j = 1, dim do inits[j] = inits[j] + (h/840) * (41 * fnc[01][j] + 216 * fnc[06][j] + 272 * fnc[07][j] + 27 * fnc[08][j] + 27 * fnc[09][j] + 36 * fnc[10][j] + 180 * fnc[11][j] + 41 * fnc[12][j]) end 

                t0 = t0 + h
             end
             return t0, inits  -- 更新した 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 = rkShanks8({xDot, yDot}, 0, {0, 0}, 0.25, _) -- (函数リスト, t0, 初期値リスト, 刻み幅 [, 内部分割数])

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

end