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

参考: Excelで学ぶ微分方程式, pp.171-176

Pend = class()
function Pend:init(funcs, t0, inits, h, numOfDiv, len, radius)
   self.len    = len
   self.radius = radius
 
   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:drawArc(sx - radiusOnScr, sy - radiusOnScr, diaOnScr, diaOnScr, 0, 360)
   gc:drawLine(cx, cy, sx, sy)
   return self
end


--- 確かめる。
do

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/2
   UNIT      = H/2.7
   FONT_SIZE = H/16
end

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

local pend
local function reset()
         pend = Pend({radDot, velDot}, 0, {math.pi, -5}, 1/64, _, 1, 0.15) -- (函数リスト, t0, {初期角度,初期角速度}, 刻み幅 [, 内部分割数], 糸の長さ, 分銅の半径)
      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:updateRK() -- タイマーの tick するたびにインスタンスを更新する。
       --:print()
   platform.window:invalidate()
end
function on.paint(gc)
   pend:drawTime(gc, 10, 0, FONT_SIZE)
       :drawPend(gc, CX, CY, UNIT)
end

end