TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 聯成振動 2 / Nyström 法 / ゴムひもの伸び縮みを表現する

CoupledVib = class()
function CoupledVib:init(funcs, t0, inits, initsDot, h, numOfDiv)
   self.fillLine = (function()
                      local POW  = math.pow
                      local ATAN = math.atan2
                      local MIN  = math.min
                      local SQRT = math.sqrt
                      local SIN  = math.sin
                      local COS  = math.cos
                      return function(gc, x1, y1, x2, y2, area) -- (gc, 始点, 終点, 面積)
                                 local angle = ATAN(y2 - y1, x2 - x1)
                                 local r     = MIN(30000, area/(2 * SQRT(POW((x2 - x1), 2) + POW((y2 - y1), 2))))
                                 local rSin = r * SIN(angle)
                                 local rCos = r * COS(angle)
                                 local p1x, p1y = x1 - rSin, y1 + rCos
                                 local p2x, p2y = x2 - rSin, y2 + rCos
                                 local p3x, p3y = x2 + rSin, y2 - rCos
                                 local p4x, p4y = x1 + rSin, y1 - rCos
                                 gc:fillPolygon({p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y})
                              end
                   end)()

   self.Unpack = unpack or table.unpack  

   self.funcs    = funcs
   self.t0       = t0
   self.inits    = inits
   self.initsDot = initsDot

   self.dim      = #funcs
   self.numOfDiv = numOfDiv or 1
   self.h        = h / self.numOfDiv

   self.f   = {{}, {}, {}, {}}
   self.tmp = {{}, {}, {}}
end
function CoupledVib:updateNY()
   for i = 1, self.numOfDiv do   
      for j = 1, self.dim do self.f[1][j] = self.funcs[j](self.Unpack(self.inits)) ; self.tmp[1][j] = self.inits[j] + self.h*2/5 * self.initsDot[j] + (2/25) * self.h * self.h *  self.f[1][j]                 end
      for j = 1, self.dim do self.f[2][j] = self.funcs[j](self.Unpack(self.tmp[1])); self.tmp[2][j] = self.inits[j] + self.h*2/3 * self.initsDot[j] + (2/9)  * self.h * self.h *  self.f[1][j]                 end
      for j = 1, self.dim do self.f[3][j] = self.funcs[j](self.Unpack(self.tmp[2])); self.tmp[3][j] = self.inits[j] + self.h*4/5 * self.initsDot[j] + (4/25) * self.h * self.h * (self.f[1][j] + self.f[2][j]) end
      for j = 1, self.dim do self.f[4][j] = self.funcs[j](self.Unpack(self.tmp[3])) end

      for j = 1, self.dim do    self.inits[j] = self.inits[j] + self.h * (self.initsDot[j] + (self.h/192) * (23 * self.f[1][j] + 75  * self.f[2][j] - 27 * self.f[3][j] + 25  * self.f[4][j])) end
      for j = 1, self.dim do self.initsDot[j] =                           self.initsDot[j] + (self.h/192) * (23 * self.f[1][j] + 125 * self.f[2][j] - 81 * self.f[3][j] + 125 * self.f[4][j])  end

      self.t0 = self.t0 + self.h
   end
   return self
end
function CoupledVib:print()
   print(self.t0, table.concat(self.inits, "\t"), table.concat(self.initsDot, "\t"))
   return self
end
function CoupledVib:drawTime(gc, x, y, fontSize, color)
   local x, y     = x or 10, y or 10
   local fontSize = fontSize or 10 
   gc:setColorRGB(color or 0x000000)
   gc:setFont("sansserif", "r", fontSize or 10)
   gc:drawString(string.format("t = %.1f", self.t0), x, y)
   return self
end
function CoupledVib:drawMassPoints(gc, x1, x2, cy, unit, w, h)
   local sx1    = x1 + unit * (self.inits[1])
   local sx2    = x2 + unit * (self.inits[2])
   local RADIUS = h / 15
   local DIA    = RADIUS + RADIUS 
   local AREA   = (1.5 * DIA)^2
   gc:setColorRGB(0x7D9EC0) -- 平衡点に縦線を引く。
   gc:setPen("thin", "dashed")
   gc:drawLine(x1, 0, x1, h)
   gc:drawLine(x2, 0, x2, h)
   gc:setPen("medium") -- バネを描く。
   gc:setColorRGB(0x1874CD) ; self.fillLine(gc, 0  , cy, sx1, cy, AREA)
   gc:setColorRGB(0x63B8FF) ; self.fillLine(gc, sx1, cy, sx2, cy, AREA)
   gc:setColorRGB(0x00868B) ; self.fillLine(gc, sx2, cy, w  , cy, AREA)
   gc:setPen("thick") -- 質点を描く。
   gc:setColorRGB(0x388E8E) ; gc:drawArc(sx1 - RADIUS, cy - RADIUS, DIA, DIA, 0, 360)
   gc:setColorRGB(0x71C671) ; gc:drawArc(sx2 - RADIUS, cy - RADIUS, DIA, DIA, 0, 360)
   return self
end


--- 確かめる。
do
platform.window:setBackgroundColor(0xE7E9DE)

local isTicking = false

local W, H, X1, X2, CY, UNIT, FONT_SIZE
function on.resize()
   W, H      = platform.window:width(), platform.window:height()
   X1        = W/3     -- m1 の平衡点
   X2        = X1 + X1 -- m2 の平衡点
   CY        = H/2
   UNIT      = W/10
   FONT_SIZE = H/16
end

local cv
local function reset()
         local m1, m2 = 2, 2 -- 質量
         local k1, k2, k3 = 10, 10, 10 -- バネ係数。値の小さいほどやわらかい。値の大きいほど固い。
         local function x1DotDot(x1, x2) return -(k1/m1) * x1 + (k2/m1) * (x2 - x1) end
         local function x2DotDot(x1, x2) return -(k2/m2) * (x2 - x1) - (k3/m2) * x2 end
         local inits    = {3, 2}
         local initsDot = {0, 0}
         cv = CoupledVib({x1DotDot, x2DotDot}, 0, inits, initsDot, 1/64, _) -- インスタンスを 1 個作る。
      end

function on.construction()
   reset()
end
function on.enterKey()
   if     isTicking == false then timer.start(cv.h * cv.numOfDiv); isTicking = true
   elseif isTicking == true  then timer.stop()                   ; isTicking = false
   end
end
function on.escapeKey()
   reset()
end
function on.timer()
   cv:updateNY()
   platform.window:invalidate()
end
function on.paint(gc)
   cv:drawMassPoints(gc, X1, X2, CY, UNIT, W, H)
     :drawTime(gc, 10, 0, FONT_SIZE, 0x000000)
end

end