参考: Excelで学ぶ微分方程式, pp.171-176
Pend = class() function Pend:init(funcs, t0, inits, h, numOfDiv, len, radius) self.len = len self.radius = radius self.Unpack = unpack or table.unpack self.Sin = math.sin self.Cos = math.cos self.funcs = funcs self.t0 = t0 self.inits = inits self.dim = #funcs self.numOfDiv = numOfDiv or 1 self.h = h / self.numOfDiv self.f = {{}, {}, {}, {}} self.temp = {{}, {}, {}} end function Pend:updateRK() for i = 1, self.numOfDiv do for j = 1, self.dim do self.f[1][j] = self.funcs[j](self.t0 , self.Unpack(self.inits)) ; self.temp[1][j] = self.inits[j] + self.h/2 * self.f[1][j] end for j = 1, self.dim do self.f[2][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[1])) ; self.temp[2][j] = self.inits[j] + self.h/2 * self.f[2][j] end for j = 1, self.dim do self.f[3][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[2])) ; self.temp[3][j] = self.inits[j] + self.h * self.f[3][j] end for j = 1, self.dim do self.f[4][j] = self.funcs[j](self.t0 + self.h , self.Unpack(self.temp[3])) end for j = 1, self.dim do self.inits[j] = self.inits[j] + self.h/6 * (self.f[1][j] + 2 * (self.f[2][j] + self.f[3][j]) + self.f[4][j]) end self.t0 = self.t0 + self.h end return self end function Pend:print() print(self.t0, "{"..table.concat(self.inits, ", ").."}") return self end function Pend: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 Pend:drawPend(gc, cx, cy, unit) local lenOnScr = unit * self.len local radiusOnScr = unit * self.radius local diaOnScr = 2 * radiusOnScr local sx = cx + lenOnScr * self.Sin(self.inits[1]) local sy = cy + lenOnScr * self.Cos(self.inits[1]) gc:drawArc(sx - radiusOnScr, sy - radiusOnScr, diaOnScr, diaOnScr, 0, 360) gc:drawLine(cx, cy, sx, sy) return self end --- 確かめる。 do local isTicking = false local W, H, CX, CY, UNIT, FONT_SIZE function on.resize() W, H = platform.window:width(), platform.window:height() CX, CY = W/2, H/2 UNIT = H/2.7 FONT_SIZE = H/16 end -- この聯立微分方程式を解く。 local function radDot(t, rad, vel) return vel end local function velDot(t, rad, vel) return -(9.8/1) * math.sin(rad) - 0.1 * vel end -- 0.1 * vel は諸々の減衰項。 local pend local function reset() pend = Pend({radDot, velDot}, 0, {math.pi, -5}, 1/64, _, 1, 0.15) -- (函数リスト, t0, {初期角度,初期角速度}, 刻み幅 [, 内部分割数], 糸の長さ, 分銅の半径) 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:updateRK() -- タイマーの tick するたびにインスタンスを更新する。 --:print() platform.window:invalidate() end function on.paint(gc) pend:drawTime(gc, 10, 0, FONT_SIZE) :drawPend(gc, CX, CY, UNIT) end end