TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 4 of 4 / 確かめる

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

function cowell4(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 staticVals(funcs, inits, h)
            local dim  = #funcs
            local w    = {}
            local f    = {{},{},{},{}}
            local temp = {}
            for i = 1, 4 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[4][j] - inits[3][j] + h*h * (204 * f[4][j] + 14 * f[3][j] + 4 * f[2][j] - f[1][j]) / 240
               temp[j] = 2 * f[4][j] - f[3][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 * 19 * temp[j]/240 -- 第 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 * 19 * f[j]/240 -- 第 2 近似値
            end
            return approx1, approx2, maxOfErr(approx1, approx2) -- {第 1 近似値}, {第 2 近似値}, 最大誤差  
         end
   local function approx2(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 * 19 * f[j]/240 -- 次の近似値
            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 = approx2(funcs, w, approx2nd, h)        -- 差がなくなるまで近似値の計算を繰り返して、差がなくなったらループを抜けて初期値を更新する。
                  count = count + 1
              end
              table.remove(inits, 1)
              inits[4] = approx2nd
              t0 = t0 + h
              return t0, inits, count -- 独立変数の値, {4 組の初期値}, 近似回数
           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

-- 4 組の初期値にはこの値を使う。 {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}
local inits = {{1.0509143, -0.4038388},
               {1.0989720,  0.0      },
               {1.0509143,  0.4038388},
               {0.9168453,  0.7754147}}
local h  =  20
local t0 = -20 

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

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

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

end

f:id:ti-nspire:20170611163708p:plain
t = 1600 のときの真値は x = -0.5893305... , y = 1.9558757...