TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 円制限三体問題 3 of 3 / トロヤ群小惑星の秤動 / 補外法

参考:


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