参考: パソコンで見る天体の動き, p.166
簡単にするため、質点同士の質量差が極端であるものとし、極端に重いほうの質点 (たとえば太陽とする) が原点から動かないものとする。長半径 1、離心率 0.999、周期 2 の条件で極端に軽いほうの質点 (何らかの天体とする) を遠日点から抛ってみる。本来は楕圓軌道を描き続けるはずであるが、近日点附近での計算誤差が大きすぎるため楕圓軌道を外れてしまう。
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("x", self.inits[1]) var.store("y", self.inits[3]) var.store("t", self.t0) var.store("div", self.numOfDiv) return self end --- 確かめる。 do -- この聯立微分方程式を解く。 local tb local function reset() local function xDot(t, x, vx, y, vy) return vx end local function vxDot(t, x, vx, y, vy) return -x / (math.sqrt(x^2 + y^2))^3 end local function yDot(t, x, vx, y, vy) return vy end local function vyDot(t, x, vx, y, vy) return -y / (math.sqrt(x^2 + y^2))^3 end local inits = { -1.999, -- x 0, -- vx 0, -- y -0.02236627199, -- vy } tb = TwoBodies({xDot, vxDot, yDot, vyDot}, 0, inits, 1/(2^6), 1e-11) 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() :print() end end
