TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 5 of 5 / 補外法 / 質点同士の接近時のみ正則化する / それに意味のないとき

参考: パソコンで見る天体の動き, p.166-175

原点に固定してある超重量物との距離が 0.65 よりも近くなったときだけ假想時間軸に移る。


TwoBodies = class()
function TwoBodies:init(funcs, t0, inits, h, numOfDiv, inWhichDomain)
   self.unpack = unpack or table.unpack
   self.Sqrt   = math.sqrt

   self.inWhichDomain = inWhichDomain   
   self.numOfDiv      = numOfDiv or 1
   self.funcs         = funcs
   self.step          = h
   self.h             = self.step / self.numOfDiv
   self.t0            = t0
   self.inits         = inits
   
   self.convert2Reals = function(r, inits)
                           local u, uPrm, v, vPrm = inits[1], inits[2], inits[3], inits[4]
                           local x = u^2 - v^2
                           local y = 2 * u * v
                           local xDot = (2/r) * (u * uPrm - v * vPrm)  
                           local yDot = (2/r) * (v * uPrm + u * vPrm)
                           return x, xDot, y, yDot
                        end
   self.setInits = function()
                      if self.inWhichDomain == "Real" then
                         self.x, self.xDot, self.y, self.yDot = unpack(self.inits)
                         self.tVirtual  = 0
                         self.tReal     = self.t0
                         self.r         = self.Sqrt((self.x - 0)^2 + (self.y - 0)^2)
                         self.initsReal = self.inits 
                      elseif self.inWhichDomain == "Virtual" then
                         self.u, self.uPrm, self.v, self.vPrm, self.tReal = unpack(self.inits)
                         self.tVirtual = self.t0
                         self.r        = self.u^2 + self.v^2
                         self.x, self.xDot, self.y, self.yDot = self.convert2Reals(self.r, {self.u, self.uPrm, self.v, self.vPrm})
                      end
                   end
   self.list2oneColMat = function(list) 
                            local oneColMat = {}
                            for i = 1, #list do
                               oneColMat[i] = {}
                               oneColMat[i][1] = list[i]
                            end
                            return oneColMat
                         end
   self.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
                       local temp = {}
                       local t    = {}
                       local x    = self.newMat(dim, 1, _)
                       local X    = self.newMat(dim, 1, _)
              
                       for s = 1, S do
                          local N  = 2^s
                          local hs = h/N
                 
                          t[1] = t0 + hs
                          for i = 1, dim
                             do x[i][1] = inits[i] + hs * funcs[i](t0, self.unpack(inits))
                          end
               
                          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], self.unpack(temp))
                             end
                          end
               
                          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
                    end
   self.neville = function(listOfLists)
                     local numOfLists = #listOfLists
                     local numOfHs    = #listOfLists[1]
                     local X          = {}
                     local extrapo    = {}
                     for n = 1, numOfLists do
                        X[n] = self.list2oneColMat(listOfLists[n])
                        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)
                           end
                        end
                        extrapo[n] = X[n][numOfHs][numOfHs]
                     end
                     return extrapo
                  end
    self.setInits()
end
function TwoBodies: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
      self.setInits()
   end
   return self
end
function TwoBodies:print()
   print(self.t0, table.concat(self.inits, "\t, "))
   return self
end


--- 確かめる。
do

local unpack = unpack or table.unpack
local Sqrt   = math.sqrt
local Ceil   = math.ceil
local Ln    = math.log

local function sign(n)
         return (n >= 0) and 1 or -1
      end
local function convert2Virtuals(t, r, inits)
         local x, xDot, y, yDot = inits[1], inits[2], inits[3], inits[4]
         local u    = Sqrt((r + x)/2)
         local v    = Sqrt((r - x)/2) * sign(y)
         local uPrm = 0.5 * ( u * xDot + v * yDot)
         local vPrm = 0.5 * (-v * xDot + u * yDot)
         return u, uPrm, v, vPrm, t
      end
local function varStore(t_real, t_virtual, x, y, r)
          var.store("t_real", t_real)
          var.store("t_virtual", t_virtual)
          var.store("x", x)
          var.store("y", y)
          var.store("r", r)
       end
local function ceilB(num)
          if (num > 0) and (num < 1) then 
             return 2^Ceil(Ln(num)/Ln(2))
          else 
             return Ceil(num)
          end
       end       

local e = 1 - (1e-14)                                                    -- 離心率
local m = {x = -(1 + e), y = 0, xDot = 0, yDot = -Sqrt(1 - e^2)/(1 + e)} -- 初期条件 
local r = Sqrt((m.x - 0)^2 + (m.y - 0)^2)                                -- 初期条件での 2 体間の距離
local E = (m.xDot^2 + m.yDot^2)/2 - (1/r)                                -- 全エネルギー
local R = 0.65                                                           -- 2 体間の距離の閾値  

local hReal    = 0.0625 
local hVirtual = ceilB(hReal/R)   
      
-- 現実時間軸で解く微分方程式
local function    xDot(t, x, xDot, y, yDot) return xDot                     end
local function xDotDot(t, x, xDot, y, yDot) return -x / (Sqrt(x^2 + y^2))^3 end
local function    yDot(t, x, xDot, y, yDot) return yDot                     end
local function yDotDot(t, x, xDot, y, yDot) return -y / (Sqrt(x^2 + y^2))^3 end
local funcsReal = {xDot, xDotDot, yDot, yDotDot}


-- 假想時間軸で解く微分方程式
local function    uPrm(t, u, uPrm, v, vPrm) return uPrm      end
local function uPrmPrm(t, u, uPrm, v, vPrm) return (E/2) * u end
local function    vPrm(t, u, uPrm, v, vPrm) return vPrm      end
local function vPrmPrm(t, u, uPrm, v, vPrm) return (E/2) * v end
local function    tPrm(t, u, uPrm, v, vPrm) return u^2 + v^2 end
local funcsVirtual = {uPrm, uPrmPrm, vPrm, vPrmPrm, tPrm}         


local tb
local function updateScreen()
         varStore(tb.tReal, tb.tVirtual, tb.x, tb.y, tb.r)
      end
local function reset()
         local initsReal    = {m.x, m.xDot, m.y, m.yDot}
         local initsVirtual = {convert2Virtuals(0, r, initsReal)}
         if     r >= R then tb = TwoBodies(funcsReal   , 0, initsReal   , hReal   , _, "Real"   )      
         elseif r <  R then tb = TwoBodies(funcsVirtual, 0, initsVirtual, hVirtual, _, "Virtual")
         end
         updateScreen()
      end

local reset_off_on = nil
function on.construction()
   reset()
   var.store("reset_off_on", 0)
   var.monitor("reset_off_on")
end
function on.varChange()
   reset_off_on = var.recall("reset_off_on")
   if     reset_off_on >  0.7 then timer.start(tb.step)
   elseif reset_off_on < -0.7 then reset()
   else                            timer.stop()
   end
end
function on.timer()
   tb:updateEX()                                                            -- 更新するたびに
   
   if (tb.inWhichDomain == "Real") and (tb.r < R) then                      -- 今現実時間軸にいて、かつ 2 体間の距離が近づきすぎたら、
      local initsVirtual = {convert2Virtuals(tb.tReal, tb.r, tb.initsReal)} -- 現実値を假想値に変換して、
      tb = TwoBodies(funcsVirtual,0 ,initsVirtual , hVirtual, _, "Virtual")  -- 假想時間軸で 2 体をインスタンス化し直すが、
   elseif (tb.inWhichDomain == "Virtual") and (tb.r >= R) then              -- 今假想時間軸にいて、かつ 2 体間の距離が十分に離れたら、
      tb = TwoBodies(funcsReal, tb.tReal, {tb.x, tb.xDot, tb.y, tb.yDot}, hReal, _, "Real"   )   -- 現実時間軸で 2 体をインスタンス化し直す。
   end 
   updateScreen()
end

end

.tns
twoBodies10.tns - Google ドライブ