TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 2 / Fehlberg 法 / それに意味のないとき / 質点同士が極端に近づくとき

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


簡単にするため、質点同士の質量差が極端であるものとし、極端に重いほうの質点 (たとえば太陽とする) が原点から動かないものとする。長半径 1、離心率 0.999、周期 2 \pi の条件で極端に軽いほうの質点 (何らかの天体とする) を遠日点から抛ってみる。本来は楕圓軌道を描き続けるはずであるが、近日点附近での計算誤差が大きすぎるため楕圓軌道を外れてしまう。

TwoBodies = class()
function TwoBodies:init(funcs, t0, inits, h, tol)
   self.step = h
   
   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
                           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 TwoBodies:updateFE()
   self.numOfDiv   = 1
   local t, a4, a5 = self.preCalc(self.numOfDiv)
   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 TwoBodies:print()
   print(self.t0, table.concat(self.inits, "\t"), self.numOfDiv)
   return self
end
function TwoBodies:varStore()
   var.store("x", self.inits[1])
   var.store("y", self.inits[3])
   var.store("t", self.t0)
   var.store("div", self.numOfDiv)
   return self
end



--- 確かめる。
do

-- この聯立微分方程式を解く。
local tb
local function reset()
         local function  xDot(t, x, vx, y, vy) return vx end
         local function vxDot(t, x, vx, y, vy) return -x / (math.sqrt(x^2 + y^2))^3 end
         local function  yDot(t, x, vx, y, vy) return vy end
         local function vyDot(t, x, vx, y, vy) return -y / (math.sqrt(x^2 + y^2))^3 end         

         local inits = {
                        -1.999, -- x
                        0, -- vx
                        0, -- y
                        -0.02236627199, -- vy
                       }
         tb = TwoBodies({xDot, vxDot, yDot, vyDot}, 0, inits, 1/(2^6), 1e-11)
         tb:varStore()
      end

local reset_off_on = nil
function on.construction()
   reset()
   var.store("reset_off_on", 0)
   var.monitor("reset_off_on")
end
function on.escapeKey()
   reset()
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:updateFE()
     :varStore()
     :print()
end

end

.tns
twoBodies2.tns - Google ドライブ