参考:
- パソコンで見る天体の動き, pp.162-166
- Excelで操る! ここまでできる科学技術計算, pp.80-86
L4 を原点に置き、太陽と木星とを結ぶ直線を x 軸に平行に固定して (つまり廻轉座標系に) 描く。ここでは両テキストと同じく Fehlberg 法を使ってみる。L4 から y 方向へ 0.01 AU ずれたところから 0.5 年刻みで 150 年間の計算する。
(全体図)
(L4 近傍の拡大図)
Fehlberg = class() function Fehlberg:init(funcs, t0, inits, h, tol) 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 -- 内部分割数が 1 のときだけ 4 次の解を計算する。 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 Fehlberg:updateFE() self.numOfDiv = 1 local t, a4, a5 = self.preCalc(self.numOfDiv) -- とりあえず元の刻み幅 (内部分割数 1) で 1 回だけ計算し、 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 Fehlberg:print() print(self.t0, table.concat(self.inits, "\t"), self.numOfDiv) return self end -- 確かめる。 do local Sqrt = math.sqrt 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 = 5.2026 -- 太陽と木星との距離〔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.01 local etaDot0 = 0 local inits = {xi0, xiDot0, eta0, etaDot0} local h = 1/2^1 local t0 = 0 local tMax = 150 local tol = 1e-5 -- 初期値を更新しながら必要な回数だけ繰り返す。 local trojan = Fehlberg(funcs, t0, inits, h, tol) local x = {trojan.inits[1]} local y = {trojan.inits[3]} for i = 1, tMax/h do trojan:updateFE() :print() x[#x+1] = trojan.inits[1] y[#y+1] = trojan.inits[3] end var.store("x", x) var.store("y", y) end
