TI-Nspire & Lua / Fehlberg 法 1 / 4 次の公式が 5 次の公式に埋め込まれた 6 段 5 次の解法を函数化する

参考: パソコンで見る天体の動き, p.103-109

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

   local dim = #funcs

   local f1  , f2  , f3  , f4  , f5   ,f6   = {}, {}, {}, {}, {}, {}
   local tmp1, tmp2, tmp3, tmp4, tmp5       = {}, {}, {}, {}, {}

   local t0     = t0
   local inits4 = {}
   local inits  = inits

   return function()
             for j = 1, dim do f1[j] = funcs[j](t0              , unpack(inits))  ; tmp1[j] = inits[j] + h * (1/4)    * f1[j]                                                                                            end
             for j = 1, dim do f2[j] = funcs[j](t0 + h * (1/4)  , unpack(tmp1))   ; tmp2[j] = inits[j] + h * (1/32)   * (3         * f1[j] + 9    * f2[j])                                                               end
             for j = 1, dim do f3[j] = funcs[j](t0 + h * (3/8)  , unpack(tmp2))   ; tmp3[j] = inits[j] + h * (1/2197) * (1932      * f1[j] - 7200 * f2[j] + 7296        * f3[j])                                         end
             for j = 1, dim do f4[j] = funcs[j](t0 + h * (12/13), unpack(tmp3))   ; tmp4[j] = inits[j] + h            * ((439/216) * f1[j] - 8    * f2[j] + (3680/513)  * f3[j] - (845/4104)  * f4[j])                   end
             for j = 1, dim do f5[j] = funcs[j](t0 + h          , unpack(tmp4))   ; tmp5[j] = inits[j] + h            * ((-8/27)   * f1[j] + 2    * f2[j] - (3544/2565) * f3[j] + (1859/4104) * f4[j] - (11/40) * f5[j]) end
             for j = 1, dim do f6[j] = funcs[j](t0 + h * (1/2)  , unpack(tmp5))                                                                                                                                          end

             for j = 1, dim do inits4[j] = inits[j] + h * ((25/216) * f1[j] + (1408/2565)  * f3[j] + (2197/4104)   * f4[j] - (1/5)  * f5[j])                  end -- 4 次の解
             for j = 1, dim do inits [j] = inits[j] + h * ((16/135) * f1[j] + (6656/12825) * f3[j] + (28561/56430) * f4[j] - (9/50) * f5[j] + (2/55) * f6[j]) end -- 5 次の解

             t0 = t0 + h


             return t0, inits4, inits
          end
end



-- 確かめる。
do

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

-- この初期値で解く。
local t0     = 0
local x0, y0 = 0, 0
local h      = 1/(2^2)

-- 1 ステップだけ計算する函数を作って、
local fe = fehlberg({xDot, yDot}, t0, {x0, y0}, h)

-- 1 ステップだけ計算する。返値は左から独立変数の値、{4 次の解}、{5 次の解}
local t, a4, a5 = fe()
print(t, "{"..table.concat(a4, ", ").."}", "{"..table.concat(a5, ", ").."}")

end

f:id:ti-nspire:20170602181928p:plain