参考:
- パソコンで見る天体の動き, pp.160-166
- Excelによる数値解析―オイラー法のうまい使い方, pp.94-105, pp.145-158
L4 を中心に、無視できるほど軽い質点をいくつか配置する。その動きを補外法で計算する。
Extrapo = class() function Extrapo:init(funcs, t0, inits, h, numOfDiv) 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:drawTrojan(gc, cx, cy, radius, unit) local sx = cx + unit * self.inits[1] local sy = cy - unit * self.inits[3] local diameter = radius + radius gc:setColorRGB(0x7171C6) gc:setPen("thin") gc:drawArc(sx - radius, sy - radius, diameter, diameter, 0, 360) return self end -- 確かめる。 do platform.window:setBackgroundColor(0xE7E9DE) local isTicking = false local Sqrt = math.sqrt local R = 5.2026 -- 太陽と木星との距離〔AU〕 local W, H, CX, CY, UNIT, FONT_SIZE, RADIUS, DIAMETER, SunSX, SunSY, JupSX, JupSY function on.resize() W, H = platform.window:width(), platform.window:height() CX, CY = W/2, H/4 UNIT = W/35 FONT_SIZE = H/16 RADIUS = H/30 DIAMETER = RADIUS + RADIUS SunSX = CX - UNIT * R/2 -- 太陽のスクリーン座標x SunSY = CY + UNIT * Sqrt(3) * R/2 -- 太陽のスクリーン座標y JupSX = CX + UNIT * R/2 -- 木星のスクリーン座標x JupSY = SunSY -- 木星のスクリーン座標y end local reset = function() local Pi = math.pi local AU = 149597870700 -- 1 天文単位〔m〕 local METER = 1 / AU -- 1 メートル〔AU〕 local YEAR = 31557600 -- 1 年〔秒〕 local SECOND = 1 / YEAR -- 1 秒〔年〕 local Gc = 6.67408e-11 -- 万有引力定数〔m^3/(kg・s^2)〕 local G = Gc * (METER^3)/(SECOND^2) -- 万有引力定数〔AU^3/(kg・year^2)〕 local S = 1.9884e30 -- 太陽の質量〔kg〕 local J = 1.8986e27 -- 木星の質量〔kg〕 local R = R -- 太陽と木星との距離〔AU〕 local P = 2 * Pi * Sqrt(R^3 / (G * (S + J))) -- 公転周期〔年〕 local n = 2 * Pi / P -- 平均運動〔1/年〕 local function r (X, Y) return Sqrt(X^2 + Y^2) end local function rmj(X, Y) return Sqrt((R - X)^2 + Y^2) end local function X (xi) return R / 2 + xi end local function Y (eta) return Sqrt(3) * R / 2 + eta end local function xiDot (t, xi, xiDot, eta, etaDot) return xiDot end local function xiDotDot (t, xi, xiDot, eta, etaDot) return 2 * n * etaDot + n^2 * X(xi) - G * S * X(xi)/r(X(xi), Y(eta))^3 + G * J * ((R - X(xi))/rmj(X(xi), Y(eta))^3 - 1/R^2) end local function etaDot (t, xi, xiDot, eta, etaDot) return etaDot end local function etaDotDot(t, xi, xiDot, eta, etaDot) return -2 * n * xiDot + n^2 * Y(eta) - G * S * Y(eta)/r(X(xi), Y(eta))^3 + G * J * (-Y(eta)/rmj(X(xi), Y(eta))^3) end local funcs = {xiDot, xiDotDot, etaDot, etaDotDot} local xi0 = 0 local xiDot0 = 0 local eta0 = 0.1 local etaDot0 = 0 local h = 1/2^6 local t0 = 0 local trojan = {} for i = 1, 50 do trojan[i] = Extrapo(funcs, t0, {xi0, xiDot0, math.random()/10 - math.random()/20, etaDot0}, h, _) end return trojan end local trojan function on.construction() trojan = reset() end function on.enterKey() if isTicking == false then timer.start(0.01); isTicking = true elseif isTicking == true then timer.stop() ; isTicking = false end end function on.escapeKey() on.construction() end function on.timer() for i = 1, #trojan do trojan[i]:updateEX() end platform.window:invalidate() end function on.paint(gc) gc:setPen("thin") gc:setFont("sansserif", "r", FONT_SIZE) gc:setColorRGB(0x000000) gc:drawString(string.format("year: %0.1f", trojan[1].t0), 10,10) gc:setColorRGB(0xC67171) gc:drawLine(SunSX - RADIUS, SunSY - RADIUS, SunSX + RADIUS, SunSY + RADIUS) gc:drawLine(SunSX + RADIUS, SunSY - RADIUS, SunSX - RADIUS, SunSY + RADIUS) gc:drawString("Sun", SunSX - DIAMETER, SunSY, "bottom") gc:setColorRGB(0x388E8E) gc:drawLine(JupSX - RADIUS, JupSY - RADIUS, JupSX + RADIUS, JupSY + RADIUS) gc:drawLine(JupSX + RADIUS, JupSY - RADIUS, JupSX - RADIUS, JupSY + RADIUS) gc:drawString("Jupiter", JupSX - DIAMETER, JupSY, "bottom") gc:setColorRGB(0x7171C6) gc:drawLine(CX - RADIUS, CY - RADIUS, CX + RADIUS, CY + RADIUS) gc:drawLine(CX + RADIUS, CY - RADIUS, CX - RADIUS, CY + RADIUS) gc:drawString("L4", CX - RADIUS, CY - RADIUS, "bottom") for i = 1, #trojan do trojan[i]:drawTrojan(gc, CX, CY, RADIUS, UNIT) end end end