参考: パソコンで見る天体の動き, 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