TI-Nspire & Lua / 補外法 5 / 補外にかける点列を中点法で求める

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

function newMat(numRows, numCols, val)
   local mat = {}
   for i = 1, numRows do
      mat[i] = {} 
      for j = 1, numCols do
         mat[i][j] = val or nil
      end
   end
   return mat
end
function midPoints(funcs, t0, inits, h, s)
   local temp = {}
   local t    = {}
   local x    = newMat(#funcs, 1, _) -- 中点法で求めた値を入れておくための 1 列行列を用意しておく。
   local X    = newMat(#funcs, 1, _) -- 補外にかける値を入れておくための 1 列行列を用意しておく。
   
   for s = 1, s do
      local N  = 2^s
      local hs = h/N -- 内部分割する刻み幅の設定はテキストどおり h/2^s とする。
      
      t[1] = t0 + hs
      for i = 1, #funcs
         do x[i][1] = inits[i] + hs * funcs[i](t0, unpack(inits))
      end -- これで x1, y1, ... が得られた (このスクリプトでの表現は x[1][1], x[2][1], ...)。
      
      for n = 0, N - 1 do
         for i = 1, #funcs do
            temp[i] = x[i][n+1]
         end
         t[n+2] = (t[n] or t0) + 2 * hs
         for i = 1, #funcs do
            x[i][n+2] = (x[i][n] or inits[i]) + 2 * hs * funcs[i](t[n+1], unpack(temp))
         end
      end -- これで {x1, x2, ..., xN+1}, {y1, y2, ..., yN+1}, ... が得られた (このスクリプトでの表現は {x[1][1], x[1][2], ..., x[1][N+1]}, {x[2][1], x[2][2], ..., x[2][N+1]})。 
      
      for i = 1, #funcs do
         X[i][s] = (1/4) * (x[i][N-1] + 2 * x[i][N] + x[i][N+1]) -- 補外にかける点列を計算する。
      end
   end
   return X -- { {X1, X2, ...}, {Y1, Y2, ...}, ... } という形で返す (このスクリプトでの表現は { {X[1][1], X[1][2], ...}, {X[2][1], X[2][2], ...}, ... })。
end



-- 確かめる。
do

-- この聯立微分方程式から、補外を適用する点列を求める。
local function xDot(t, x, y) return y     end     
local function yDot(t, x, y) return t - x end

local a = midPoints({xDot, yDot}, 0, {0, 0}, 0.5, 7)
for i = 1, #a do
   print("{"..table.concat(a[i], ", ").."}")
end

end

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