TI-Nspire & Lua / 特殊な方程式に対する Nyström 法

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

特殊な方程式に対する Nyström 法を使って二階のまま解く。ここでいう「特殊な方程式」とは、一次微分も独立変数も含まない二階の方程式のことである。

function nystroem(funcs, t0, inits, initsDot, h, numOfDiv)
   local unpack = unpack or table.unpack  

   local dim      = #funcs
   local numOfDiv = numOfDiv or 1
   local h        = h / numOfDiv

   local f1  , f2  , f3  , f4   = {}, {}, {}, {}
   local tmp1, tmp2, tmp3       = {}, {}, {}

   local t0       = t0
   local inits    = inits
   local initsDot = initsDot

   return function()
             for i = 1, numOfDiv do   
                for j = 1, dim do f1[j] = funcs[j](unpack(inits)) ; tmp1[j] = inits[j] + (2/5) * h * initsDot[j] + (2/25) * h * h *  f1[j]          end
                for j = 1, dim do f2[j] = funcs[j](unpack(tmp1))  ; tmp2[j] = inits[j] + (2/3) * h * initsDot[j] + (2/9)  * h * h *  f1[j]          end
                for j = 1, dim do f3[j] = funcs[j](unpack(tmp2))  ; tmp3[j] = inits[j] + (4/5) * h * initsDot[j] + (4/25) * h * h * (f1[j] + f2[j]) end
                for j = 1, dim do f4[j] = funcs[j](unpack(tmp3))                                                                                    end

                for j = 1, dim do    inits[j] =  inits[j] + h * (initsDot[j] + (h/192) * (23 * f1[j] + 75  * f2[j] - 27 * f3[j] + 25  * f4[j])) end
                for j = 1, dim do initsDot[j] =                  initsDot[j] + (h/192) * (23 * f1[j] + 125 * f2[j] - 81 * f3[j] + 125 * f4[j])  end

                t0 = t0 + h
             end
             return t0, unpack(inits) -- 多値を返す場合の一般的な注意: unpack(inits), t0 の順番で指定すると、多値を返すはずの unpack 函数は最初の要素しか返さなくなる。
             
          end
end


-- 振子の運動方程式で確かめる。
do
local isTicking = false

-- この二階微分方程式を解く。
local function radDotDot(rad) return -9.80665 * math.sin(rad) end

-- この初期値で解く。
local t0   = 0
local rad0 = math.pi * 179.9/180
local vel0 = 0
local h    = 1/(2^6)

-- 1 ステップだけ計算する函数を作って、
local ny = nystroem({radDotDot}, t0, {rad0}, {vel0}, h, _)

-- タイマーの tick するたびにその函数を実行して初期値を更新する。
function on.timer()
   t0, rad0 = ny()
   --print(t0, rad0) -- コンソールでの確認用。
   platform.window:invalidate()
end

-- Enter キーでタイマーを toggle する。
function on.enterKey()
   if   isTicking == false then timer.start(h) isTicking = true
   else                         timer.stop()   isTicking = false
   end
end

-- 振子の棒を描く。
function on.paint(gc, _, _, W, H)
   gc:drawString("sec: "..string.format("%.1f", t0                ), 0, 0   )
   gc:drawString("deg: "..string.format("%.1f", 180 * rad0/math.pi), 0, H/10)
   gc:drawLine(W/2, H/2, W/2 + H * 0.45 * math.sin(rad0), H/2 + H * 0.45 * math.cos(rad0))
end

end