Lua の実行画面は使わずに Graphs アプリケーションに描画してみる。再生、一時停止、リセットの各操作はスライダーをボタン代わりにして行う。
TwoBodies = class() function TwoBodies:init(funcs, t0, inits, h, tol) self.step = h self.funcs = funcs self.t0 = t0 self.inits = inits self.h = h self.numOfDiv = false self.dim = #funcs self.tol = tol or 0.001 self.Unpack = unpack or table.unpack self.Ln = math.log self.Floor = math.floor self.Abs = math.abs self.Max = math.max self.floorB = function(num) if (num > 0) and (num < 1) then return 2^self.Floor(self.Ln(num)/self.Ln(2)) else return self.Floor(num) end end self.hNew = function(h, err, tol) if err > tol then return 0.9 * h * (tol * h/(h * err))^(1/4) else return h end end self.maxOfErr = function(listA, listB) local sute = {} for i = 1, #listA do sute[i] = self.Abs(listA[i] - listB[i]) end return self.Max(self.Unpack(sute)) end self.preCalc = function(numOfDiv) local f = {{}, {}, {}, {}, {}, {}} local tmp = {{}, {}, {}, {}, {}} local inits4 = {} local preInits = {self.Unpack(self.inits)} local preT0 = self.t0 local preH = self.h/numOfDiv for i = 1, numOfDiv do for j = 1, self.dim do f[1][j] = self.funcs[j](preT0 , self.Unpack(preInits)); tmp[1][j] = preInits[j] + preH * (1/4) * f[1][j] end for j = 1, self.dim do f[2][j] = self.funcs[j](preT0 + preH * (1/4) , self.Unpack(tmp[1])) ; tmp[2][j] = preInits[j] + preH * (1/32) * (3 * f[1][j] + 9 * f[2][j]) end for j = 1, self.dim do f[3][j] = self.funcs[j](preT0 + preH * (3/8) , self.Unpack(tmp[2])) ; tmp[3][j] = preInits[j] + preH * (1/2197) * (1932 * f[1][j] - 7200 * f[2][j] + 7296 * f[3][j]) end for j = 1, self.dim do f[4][j] = self.funcs[j](preT0 + preH * (12/13), self.Unpack(tmp[3])) ; tmp[4][j] = preInits[j] + preH * ((439/216) * f[1][j] - 8 * f[2][j] + (3680/513) * f[3][j] - (845/4104) * f[4][j]) end for j = 1, self.dim do f[5][j] = self.funcs[j](preT0 + preH , self.Unpack(tmp[4])) ; tmp[5][j] = preInits[j] + preH * ((-8/27) * f[1][j] + 2 * f[2][j] - (3544/2565) * f[3][j] + (1859/4104) * f[4][j] - (11/40) * f[5][j]) end for j = 1, self.dim do f[6][j] = self.funcs[j](preT0 + preH * (1/2) , self.Unpack(tmp[5])) end if numOfDiv == 1 then for j = 1, self.dim do inits4[j] = preInits[j] + preH * ((25/216) * f[1][j] + (1408/2565) * f[3][j] + (2197/4104) * f[4][j] - (1/5) * f[5][j]) end end for j = 1, self.dim do preInits[j] = preInits[j] + preH * ((16/135) * f[1][j] + (6656/12825) * f[3][j] + (28561/56430) * f[4][j] - (9/50) * f[5][j] + (2/55) * f[6][j]) end preT0 = preT0 + preH end return preT0, inits4, preInits end end function TwoBodies:updateFE() self.numOfDiv = 1 local t, a4, a5 = self.preCalc(self.numOfDiv) local err = self.maxOfErr(a4, a5) if err < self.tol then self.t0, self.inits = t, a5 else self.numOfDiv = self.h/self.floorB(self.hNew(self.h, err, self.tol)) t, a4, a5 = self.preCalc(self.numOfDiv) self.t0, self.inits = t, a5 end return self end function TwoBodies:print() print(self.t0, table.concat(self.inits, "\t"), self.numOfDiv) return self end function TwoBodies:varStore() var.store("x1", self.inits[1]) var.store("y1", self.inits[3]) var.store("x2", self.inits[5]) var.store("y2", self.inits[7]) var.store("t", self.t0) var.store("div", self.numOfDiv) return self end --- 確かめる。 do -- この聯立微分方程式を解く。 local tb local function reset() local m1 = 2 local m2 = 1 local function x1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vx1 end local function vx1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return -m2 * (x1 - x2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end local function y1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vy1 end local function vy1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return -m2 * (y1 - y2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end local function x2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vx2 end local function vx2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return m1 * (x1 - x2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end local function y2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vy2 end local function vy2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return m1 * (y1 - y2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end local inits = { -1, -- x1 0, -- vx1 0, -- y1 -0.2, -- vy1 1, -- x2 0, -- vx2 0, -- y2 0.4, --vy2 } tb = TwoBodies({x1Dot, vx1Dot, y1Dot, vy1Dot, x2Dot, vx2Dot, y2Dot, vy2Dot}, 0, inits, 1/(2^6), 1e-12) tb:varStore() end local reset_off_on = nil function on.construction() reset() var.store("reset_off_on", 0) var.monitor("reset_off_on") end function on.escapeKey() reset() end function on.varChange() reset_off_on = var.recall("reset_off_on") if reset_off_on == 1 then timer.start(tb.step) end if reset_off_on == 0 then timer.stop() end if reset_off_on == -1 then reset() end end function on.timer() tb:updateFE() :varStore() end end