TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 1 / Fehlberg 法 / それに意味のないとき

Lua の実行画面は使わずに Graphs アプリケーションに描画してみる。再生、一時停止、リセットの各操作はスライダーをボタン代わりにして行う。

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("x1", self.inits[1])
   var.store("y1", self.inits[3])
   var.store("x2", self.inits[5])
   var.store("y2", self.inits[7])
   var.store("t", self.t0)
   var.store("div", self.numOfDiv)
   return self
end



--- 確かめる。
do

-- この聯立微分方程式を解く。
local tb
local function reset()
         local m1 = 2
         local m2 = 1

         local function  x1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vx1 end
         local function vx1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return -m2 * (x1 - x2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end

         local function  y1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vy1 end     
         local function vy1Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return -m2 * (y1 - y2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end

         local function  x2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vx2 end
         local function vx2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return m1 * (x1 - x2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end

         local function  y2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return vy2 end     
         local function vy2Dot(t, x1, vx1, y1, vy1, x2, vx2, y2, vy2) return m1 * (y1 - y2) / (math.sqrt((x1 - x2)^2 + (y1 - y2)^2))^3 end

         local inits = {
                        -1, -- x1
                        0, -- vx1
                        0, -- y1
                        -0.2, -- vy1
                        1, -- x2
                        0, -- vx2
                        0, -- y2
                        0.4, --vy2
                       }
         tb = TwoBodies({x1Dot, vx1Dot, y1Dot, vy1Dot, x2Dot, vx2Dot, y2Dot, vy2Dot}, 0, inits, 1/(2^6), 1e-12)
         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()
end

end

.tns:
twoBodies1.tns - Google ドライブ