TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / Pendulum Waves / 古典的 Runge-Kutta 法

長さの違う複数の振子を同じ初期角から一斉に振る。

各糸の長さは、たとえば同じ 30 秒間に振子 1 は 15 往復、振子 2 は 16 往復、振子 3 は 17 往復 ...... するように設定する。

参考:

Pendulum Waves


Pend = class()
function Pend:init(funcs, t0, inits, h, numOfDiv, len, radius, color)
   self.len    = len
   self.radius = radius
   self.color  = color
 
   self.Unpack = unpack or table.unpack
   self.Sin    = math.sin
   self.Cos    = math.cos

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

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

   self.f    = {{}, {}, {}, {}}
   self.temp = {{}, {}, {}}
end
function Pend:updateRK()
   for i = 1, self.numOfDiv do   
      for j = 1, self.dim do self.f[1][j] = self.funcs[j](self.t0           , self.Unpack(self.inits))   ; self.temp[1][j] = self.inits[j] + self.h/2 * self.f[1][j] end
      for j = 1, self.dim do self.f[2][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[1])) ; self.temp[2][j] = self.inits[j] + self.h/2 * self.f[2][j] end
      for j = 1, self.dim do self.f[3][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[2])) ; self.temp[3][j] = self.inits[j] + self.h   * self.f[3][j] end
      for j = 1, self.dim do self.f[4][j] = self.funcs[j](self.t0 + self.h  , self.Unpack(self.temp[3]))                                                             end

      for j = 1, self.dim do self.inits[j] = self.inits[j] + self.h/6 * (self.f[1][j] + 2 * (self.f[2][j] + self.f[3][j]) + self.f[4][j]) end 

      self.t0 = self.t0 + self.h
   end
   return self
end
function Pend:print()
   print(self.t0, "{"..table.concat(self.inits, ", ").."}")
   return self
end
function Pend: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 Pend:drawPend(gc, cx, cy, unit)
   local lenOnScr    = unit * self.len
   local radiusOnScr = unit * self.radius
   local diaOnScr    = 2 * radiusOnScr
   
   local sx  = cx + lenOnScr * self.Sin(self.inits[1])
   local sy  = cy + lenOnScr * self.Cos(self.inits[1])
   
   gc:setColorRGB(self.color)
   gc:drawArc(sx - radiusOnScr, sy - radiusOnScr, diaOnScr, diaOnScr, 0, 360)
   gc:drawLine(cx, cy, sx, sy)
   return self
end


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

local isTicking = false

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

local colors   = {0xc6c6ff,0xc1c1ff,0xbcbcff,0xb7b7ff,0xb2b2ff,0xa8a8ff,0xa8a8ff,0xa3a3ff,0x9e9eff,0x9999ff,0x9393ff,0x8e8eff,0x8989ff,0x8484ff,0x7f7fff}
local lens     = {0.96285212588934,0.84625676136746,0.74945046484677,0.66851360342403,0.60015560057340,0.54160430424661,0.49120194251045,0.44751739096444,0.40962625507320,0.37611413577122,0.34662680356814,0.32056136826120,0.29711721061937,0.27636562941545,0.25760882741459}
local radiuses = {0.05,0.045,0.041,0.036,0.032,0.030,0.027,0.024,0.022,0.020,0.017,0.016,0.014,0.013,0.011}

-- この聯立微分方程式を解く。
local function radDot(t, rad, vel, len) return vel end     
local function velDot(t, rad, vel, len) return -(9.80665/len) * math.sin(rad) - 0.01 * vel end -- 0.01 * vel は諸々の減衰項。
local function lenDot(t, rad, vel, len) return 0 end

-- 長さの違う複数の振子をインスタンス化する。
local pends = {}
local function reset()
         for i = 1, #lens do
            pends[i] = Pend({radDot, velDot, lenDot}, 0, {0.5, 0, lens[i]}, 1/64, _, lens[i], radiuses[i], colors[i]) -- (函数リスト, t0, {初期角度,初期角速度}, 刻み幅 [, 内部分割数], 糸の長さ, 分銅の半径, 色)
         end
      end

function on.construction()
   reset()
end
function on.enterKey()
   if     isTicking == false then timer.start(pends[1].h) ; isTicking = true
   elseif isTicking == true  then timer.stop()            ; isTicking = false
   end
end
function on.escapeKey()
   reset()
end
function on.timer()
   for i = 1, #lens do
      pends[i]:updateRK() -- タイマーの tick するたびにインスタンスを更新する。
   end
   platform.window:invalidate()
end
function on.paint(gc)
   gc:setColorRGB(colors[1])
   pends[1]:drawTime(gc, 10, 0, FONT_SIZE)
   
   for i = #lens, 1, -1 do -- 奥の振子から順番に描く。
      pends[i]:drawPend(gc, CX, CY, UNIT)
   end
end

end