参考: ばね振り子 - epii's physics notes
Extrapo = class() function Extrapo:init(funcs, t0, inits, h, numOfDiv) self.Sin = math.sin self.Cos = math.cos local unpack = unpack or table.unpack self.numOfDiv = numOfDiv or 1 self.funcs = funcs self.step = h self.h = self.step / self.numOfDiv self.t0 = t0 self.inits = inits local list2oneColMat = function(list) local oneColMat = {} for i = 1, #list do oneColMat[i] = {} oneColMat[i][1] = list[i] end return oneColMat end local newMat = function(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 self.midPoints = function(funcs, t0, inits, h) local dim = #funcs local S = 7 -- テキストどおり 7 通り計算する。 local temp = {} local t = {} local x = newMat(dim, 1, _) -- 中点法で求めた値を入れておくための 1 列行列を用意しておく。 local X = newMat(dim, 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, dim 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, dim do temp[i] = x[i][n+1] end t[n+2] = (t[n] or t0) + 2 * hs for i = 1, dim 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, dim 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 self.neville = function(listOfLists) local numOfLists = #listOfLists local numOfHs = #listOfLists[1] local X = {} local extrapo = {} for n = 1, numOfLists do X[n] = list2oneColMat(listOfLists[n]) -- 既知の点列を 1 列行列に変換する。 for i = 2, numOfHs do for j = 2, i do X[n][i][j] = X[n][i][j-1] + (X[n][i][j-1] - X[n][i-1][j-1]) / (4^(j-1) - 1) -- 順次 1/2 刻みで計算する場合の公式。 end end extrapo[n] = X[n][numOfHs][numOfHs] -- 補外値だけ返す。 end return extrapo end end function Extrapo:updateEX() for i = 1, self.numOfDiv do self.inits = self.neville(self.midPoints(self.funcs, self.t0, self.inits, self.h)) self.t0 = self.t0 + self.h end return self end function Extrapo:print() print(self.t0, table.concat(self.inits, "\t")) return self end function Extrapo:drawTime(gc, x, y, fontSize) local x, y = x or 10, y or 10 local fontSize = fontSize or 10 gc:setFont("sansserif", "r", fontSize) gc:drawString(string.format("t = %.1f", self.t0), x, y) return self end function Extrapo:drawPend(gc, cx, cy, unit, pendRadius) platform.window:setBackgroundColor(0xE7E9DE) local pendDia = pendRadius * 2 local rad = self.inits[1] local slen = unit * self.inits[3] local sx = cx + slen * self.Sin(rad) local sy = cy + slen * self.Cos(rad) gc:setPen("medium") gc:setColorRGB(0x1874CD) gc:drawLine(cx, cy, sx, sy) gc:setPen("thick") gc:setColorRGB(0x71C671) gc:drawArc(sx - pendRadius, sy - pendRadius, pendDia, pendDia, 0, 360) return self end -- 確かめる。 do local isTicking = false local W, H, CX, CY, UNIT, FONT_SIZE, PEND_RADIUS function on.resize() W, H = platform.window:width(), platform.window:height() CX, CY = W/2, H/2 UNIT = H/12 FONT_SIZE = H/16 PEND_RADIUS = W/30 end local Sin = math.sin local Cos = math.cos local g = 9.80665 local m = 1.5 local len0 = 2 -- バネの自然長 local k = 20 -- この聯立微分方程式を補外法で解く。 local function radDot (t, rad, radDot, len, lenDot) return radDot end local function radDotDot(t, rad, radDot, len, lenDot) return (-2 * lenDot * radDot - g * Sin(rad)) / len end local function lenDot (t, rad, radDot, len, lenDot) return lenDot end local function lenDotDot(t, rad, radDot, len, lenDot) return len * radDot^2 + g * Cos(rad) - k * (len - len0) / m end local funcs = {radDot, radDotDot, lenDot, lenDotDot} -- この初期値で解く。 local radInit = 3 local radDotInit = 0 local lenInit = len0 * 0.1 -- バネを圧縮した状態から始める。 local lenDotInit = 0 local Inits = {radInit, radDotInit, lenInit, lenDotInit} local pend local function reset() pend = Extrapo(funcs, 0, Inits, 1/2^6, _) end function on.construction() reset() end function on.enterKey() if isTicking == false then timer.start(pend.h) ; isTicking = true elseif isTicking == true then timer.stop() ; isTicking = false end end function on.escapeKey() reset() end function on.timer() pend:updateEX() platform.window:invalidate() end function on.paint(gc) pend:drawTime(gc, 10, 0, FONT_SIZE) :drawPend(gc, CX, CY, UNIT, PEND_RADIUS) end end