TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / ばね振り子 2 of 2 / 軌跡を描く / 補外法

Extrapo = class()
function Extrapo:init(funcs, t0, inits, h, numOfDiv)

   self.Sin = math.sin
   self.Cos = math.cos
   
   self.historicData = {}
               
   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
   
   self.storeData =
   function(x, y, list, numOfPairs)
      list[#list+1] = x
      list[#list+1] = y
      if #list > numOfPairs * 2 then
         table.remove(list, 1)
         table.remove(list, 1)
      end
   end
   
   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:drawTime(gc, x, y, fontSize)
   local x, y     = x or 10, y or 10
   local fontSize = fontSize or 10 
   gc:setFont("sansserif", "r", fontSize)
   gc:drawString(string.format("t = %.1f", self.t0), x, y)
   return self
end
function Extrapo:drawPend(gc, cx, cy, unit, pendRadius)
   platform.window:setBackgroundColor(0xE7E9DE)
   local pendDia = pendRadius * 2
   self.rad     = self.inits[1]
   self.slen    = unit * self.inits[3]
   self.sx      = cx + self.slen * self.Sin(self.rad)
   self.sy      = cy + self.slen * self.Cos(self.rad)
   gc:setPen("medium")
   gc:setColorRGB(0x1874CD)
   gc:drawLine(cx, cy, self.sx, self.sy)
   gc:setPen("thick")
   gc:setColorRGB(0x71C671)
   gc:drawArc(self.sx - pendRadius, self.sy - pendRadius, pendDia, pendDia, 0, 360)
   return self
end
function Extrapo:drawLocus(gc, numOfPairs, color)
   self.storeData(self.sx, self.sy, self.historicData, numOfPairs)
   gc:setColorRGB(color or 0x000000)
   gc:setPen("thin")
   gc:drawPolyLine(self.historicData)
   return self
end


-- 確かめる。
do
local isTicking = false

local W, H, CX, CY, UNIT, FONT_SIZE, PEND_RADIUS 
function on.resize()
   W, H        = platform.window:width(), platform.window:height()
   CX, CY      = W/2, H/3
   UNIT        = H/16
   FONT_SIZE   = H/16
   PEND_RADIUS = W/30
end

local Sin = math.sin
local Cos = math.cos
local g    = 9.80665
local m    = 1.5
local len0 = 3 -- バネの自然長
local k    = 10

-- この聯立微分方程式を補外法で解く。
local function radDot   (t, rad, radDot, len, lenDot) return radDot end
local function radDotDot(t, rad, radDot, len, lenDot) return (-2 * lenDot * radDot - g * Sin(rad)) / len end
local function lenDot   (t, rad, radDot, len, lenDot) return lenDot end
local function lenDotDot(t, rad, radDot, len, lenDot) return len * radDot^2 + g * Cos(rad) - k * (len - len0) / m end
local funcs = {radDot, radDotDot, lenDot, lenDotDot}

-- この初期値で解く。
local radInit    = 3.1
local radDotInit = 0
local lenInit    = len0 * 0.2
local lenDotInit = 0
local Inits = {radInit, radDotInit, lenInit, lenDotInit}

local pend
local function reset()
         pend = Extrapo(funcs, 0, Inits, 1/2^6, _)
      end

function on.construction()
   reset()
end
function on.enterKey()
   if     isTicking == false then timer.start(pend.h) ; isTicking = true
   elseif isTicking == true  then timer.stop()        ; isTicking = false
   end
end
function on.escapeKey()
   reset()
end
function on.timer()
   pend:updateEX()
   platform.window:invalidate()
end
function on.paint(gc)
   pend:drawTime(gc, 10, 0, FONT_SIZE)
       :drawPend(gc, CX, CY, UNIT, PEND_RADIUS)
       :drawLocus(gc, 500, 0x555555)
end

end