TI-Nspire & Lua / 特殊な方程式に対する 7 段の Cowell 法

参考: パソコンで見る天体の動き, pp.112-119

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

   local function maxOfErr(listA, listB)
            local sute = {}
            for i = 1, #listA do
               sute[i] = math.abs(listA[i] - listB[i])
            end
            return math.max(unpack(sute))
         end
   local function maxOfErr(listA, listB)
            local sute = {}
            for i = 1, #listA do
               sute[i] = math.abs(listA[i] - listB[i])
            end
            return math.max(unpack(sute))
         end
   local function staticVals(funcs, inits, h)
            local dim  = #funcs
            local w    = {}
            local f    = {{},{},{},{},{},{},{}}
            local temp = {}
            for i = 1, 7 do
               for j = 1, dim do
                  f[i][j] = funcs[j](unpack(inits[i]))
               end
            end  
            for j = 1, dim do
               w[j]    = 2 * inits[7][j] - inits[6][j] + h*h * (55324 * f[7][j] - 6297 * f[6][j] + 14598 * f[5][j] - 11477 * f[4][j] + 5568 * f[3][j] - 1551 * f[2][j] + 190 * f[1][j]) / 60480
               temp[j] = 2 * f[7][j] - f[6][j]
            end
            return w, temp -- {変化しない部分}, {假定値}
         end
   local function approx(funcs, w, temp, h) -- ({微分方程式}, {変化しない部分}, {假定値}, 刻み幅)
            local dim     = #funcs
            local approx1 = {}
            local approx2 = {}
            local f       = {}
            for j = 1, dim do
               approx1[j] = w[j] + h*h * 4125 * temp[j]/60480 -- 第 1 近似値
            end
            for j = 1, dim do
               f[j] = funcs[j](unpack(approx1)) 
            end
            for j = 1, dim do
               approx2[j] = w[j] + h*h * 4125 * f[j]/60480 -- 第 2 近似値
            end
            return approx1, approx2, maxOfErr(approx1, approx2) -- {第 1 近似値}, {第 2 近似値}, 最大誤差  
         end
   local function APPROX(funcs, w, preApprox, h) -- ({微分方程式}, {変化しない部分}, {前の近似値}, 刻み幅)
            local dim    = #funcs
            local f      = {}
            local approx = {}
            for j = 1, dim do
               f[j] = funcs[j](unpack(preApprox)) 
            end
            for j = 1, dim do
               approx[j] = w[j] + h*h * 4125 * f[j]/60480 -- 次の近似値
            end
            return approx, maxOfErr(preApprox, approx) -- {次の近似値}, (前の近似値と次の近似値との最大誤差)  
         end

         return function()
                   local count = 0 
                   local w, temp = staticVals(funcs, inits, h)                 -- まず変化しない部分と最初の假定値とを求めて、
                   local approx1st, approx2nd, err = approx(funcs, w, temp, h) -- 最初の第 1 近似値、第 2 近似値、両者の差を求めて、
                   while(err > 1e-15) do                                       -- 差が存在しなかったらループを飛ばして初期値を更新するが、差が存在していたらループに入って、
                      approx2nd, err = APPROX(funcs, w, approx2nd, h)          -- 差がなくなるまで近似値の計算を繰り返して、差がなくなったらループを抜けて初期値を更新する。
                      count = count + 1
                   end
                   table.remove(inits, 1)
                   inits[7] = approx2nd
                   t0 = t0 + h
                   return t0, inits, count -- 独立変数の値, {7 組の初期値}, 近似回数
                end
end


--------------------------------------------------------------------------------
-- 確かめる
do

-- この聯立二階微分方程式で確かめる。
local function xDotDot(x, y) return -2.95912208e-4 * x / (math.sqrt(x^2 + y^2)^3) end
local function yDotDot(x, y) return -2.95912208e-4 * y / (math.sqrt(x^2 + y^2)^3) end

-- 7 組の初期値にはこの値を使う。 {{x1, y1}, {x2, y2}, ..., {x6, y6}, {x7, y7}}
local inits = {{1.0509145, -0.4038387},
               {1.0989720,  0.0      },
               {1.0509145,  0.4038387},
               {0.9168459,  0.7754143},
               {0.7200885,  1.0952558},
               {0.4849475,  1.3582345},
               {0.2301295,  1.5678560}}
local h  =  20
local t0 = -20 

-- 1 ステップだけ計算する函数を作って、
local nextInits = cowell7({xDotDot, yDotDot}, t0, inits, h)

-- 必要な回数だけその函数を実行する。
local t, sute, count
for i = 1, 81 do 
   t, sute, count = nextInits()
end

-- 独立変数の値, 近似回数, {更新された 7 組の初期値} を print して確かめる
local str = {}
for i = 1, 7 do
   str[i] = "{"..table.concat(sute[i], ", ").."}"
end
print(t, count, "{"..table.concat(str, ", ").."}")

end

t = 1600 のときの真値は x = -0.5893305... , y = 1.9558757...