TI-Nspire & Lua / Fehlberg 法 5 / 指定許容範囲に収まるかどうかを確かめる

function fehlberg(funcs, t0, inits, h, tol)
   local unpack = unpack or table.unpack
   local Ln     = math.log
   local Floor  = math.floor
   local Abs    = math.abs
   local Max    = math.max  

   local dim = #funcs

   local t0    = t0
   local inits = inits

   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 f1  , f2  , f3  , f4  , f5   ,f6   = {}, {}, {}, {}, {}, {}
            local tmp1, tmp2, tmp3, tmp4, tmp5       = {}, {}, {}, {}, {}

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

            for i = 1, numOfDiv do
               for j = 1, dim do f1[j] = funcs[j](preT0                 , unpack(preInits)); tmp1[j] = preInits[j] + preH * (1/4)    * f1[j]                                                                                            end 
               for j = 1, dim do f2[j] = funcs[j](preT0 + preH * (1/4)  , unpack(tmp1))    ; tmp2[j] = preInits[j] + preH * (1/32)   * (3         * f1[j] + 9    * f2[j])                                                               end
               for j = 1, dim do f3[j] = funcs[j](preT0 + preH * (3/8)  , unpack(tmp2))    ; tmp3[j] = preInits[j] + preH * (1/2197) * (1932      * f1[j] - 7200 * f2[j] + 7296        * f3[j])                                         end
               for j = 1, dim do f4[j] = funcs[j](preT0 + preH * (12/13), unpack(tmp3))    ; tmp4[j] = preInits[j] + preH            * ((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](preT0 + preH          , unpack(tmp4))    ; tmp5[j] = preInits[j] + preH            * ((-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](preT0 + preH * (1/2)  , unpack(tmp5))                                                                                                                                                 end

               if numOfDiv == 1 then -- 内部分割数が 1 のときだけ 4 次の解を計算する。
                  for j = 1, dim do   inits4[j] = preInits[j] + preH * ((25/216) * f1[j] + (1408/2565)  * f3[j] + (2197/4104)   * f4[j] - (1/5)  * f5[j]) end
               end
               for j = 1, dim do preInits[j] = preInits[j] + preH * ((16/135) * f1[j] + (6656/12825) * f3[j] + (28561/56430) * f4[j] - (9/50) * f5[j] + (2/55) * f6[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
          end
   
end



-- 確かめる
do

-- この微分方程式を解く。
local function xDot(t, x) return x^2 - t^2 - 2 * t + 2 end

-- この初期値で解く。
local t0 = 0
local x0 = 0
local h  = 1

-- 指定許容値ごとに 1 ステップだけ計算する函数を作って比べる。
local fe1 = fehlberg({xDot}, t0, {x0}, h, 1e-1)
local fe2 = fehlberg({xDot}, t0, {x0}, h, 1e-5)
local fe3 = fehlberg({xDot}, t0, {x0}, h, 1e-10)

-- 1 ステップだけ計算する。
-- 独立変数の値, {5 次の解}, 1 ステップの分割数, 真値 (1.5) との差を print して確かめる。
local t, a5, div = fe1()
print(t, "{"..table.concat(a5, ", ").."}", div, (unpack(a5) - 1.5))
local t, a5, div = fe2()
print(t, "{"..table.concat(a5, ", ").."}", div, (unpack(a5) - 1.5))
local t, a5, div = fe3()
print(t, "{"..table.concat(a5, ", ").."}", div, (unpack(a5) - 1.5))


end