TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 4 / 補外法 / 正則化した方程式を解く / それに意味のないとき

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

Levi-Civita の変換で正則化した方程式を解く。假想時間軸の進みは一様であるが、現実時間軸の進みは、質点同士 (一方の質点は原点に固定してある) の近づいたときに遅くなる。誤差を見積もらない限り微分方程式が解けたことにはならないが、離心率が極端に 1 に近い楕圓軌道でも計算は破綻しない。

長半径 1 の楕圓軌道を周期 2 Pi で周回するように初期条件を設定する。

TwoBodiesVirtual = class()
function TwoBodiesVirtual:init(funcs, t0, inits, h, numOfDiv)
     
   self.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
   
   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
end
function TwoBodiesVirtual: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 TwoBodiesVirtual:print()
   print(self.t0, table.concat(self.inits, "\t, "))
   return self
end


--- 確かめる。
do

local unpack = unpack or table.unpack
local Sqrt   = math.sqrt

local function sign(n)
         return (n >= 0) and 1 or -1
      end
local function distance(...)
          local aug = {...}
          if #aug == 4 then
             local x1, y1, x2, y2 = aug[1], aug[2], aug[3], aug[4]
             return Sqrt((aug[1] - aug[3])^2 + (aug[2] - aug[4])^2)
          else
             local u, v = aug[1], aug[2]
             return u^2 + v^2
          end 
       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 convert2Reals(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
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 e = 1 - (1e-14) -- 離心率
local m = {x = -(1 + e), y = 0, xDot = 0, yDot = -Sqrt(1 - e^2)/(1 + e)} -- 初期条件 
local E = (m.xDot^2 + m.yDot^2)/2 - (1/distance(m.x, m.y, 0, 0))         -- 全エネルギー      
      
--[[ 現実時間軸で解く微分方程式
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 reset()
         local initsReal    = {m.x, m.xDot, m.y, m.yDot}
         local r            = distance(initsReal[1], initsReal[3], 0, 0)
         local initsVirtual = {convert2Virtuals(0, r, initsReal)} -- 現実時間軸の初期値を假想時間軸の初期値に変換する。
         tb = TwoBodiesVirtual(funcsVirtual, 0, initsVirtual, 1/2^5, _)  -- 假想時間軸で 2 体をインスタンス化する。
         
         r = distance(tb.inits[1], tb.inits[3])
         local x, xDot, y, yDot = convert2Reals(r, tb.inits)      -- 假想値を現実値に変換する。
         varStore(tb.inits[5], tb.t0, x, y, r)
      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 ==  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:updateEX()
   local r          = distance(tb.inits[1], tb.inits[3])
   local x, _, y, _ = convert2Reals(r, tb.inits) -- 假想値を現実値に変換する。
   varStore(tb.inits[5], tb.t0, x, y, r)
end

end

.tns
twoBodies7.tns - Google ドライブ