TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 聯成振動 1 / Nyström 法 / それに意味のないとき

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

テキストは古典的 Runge-Kutta 法を使っているが、ここでは、特殊な方程式に対する Nyström 法を使ってみる。

縦の点線が平衡点:

CoupledVib = class()
function CoupledVib:init(funcs, t0, inits, initsDot, h, numOfDiv)
   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 
   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) ; gc:drawLine(0  , cy, sx1, cy)
   gc:setColorRGB(0x63B8FF) ; gc:drawLine(sx1, cy, sx2, cy)
   gc:setColorRGB(0x00868B) ; gc:drawLine(sx2, cy, w  , cy)
   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/5
   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    = {1, 0}
         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

Excelで学ぶ微分方程式

―――――――――――――――――――
厳密解が求まるので近似解を求める意味は計算の確認以外にはない。上と同じ条件で厳密解を求めると下のようになる。

DSolve[{x1''[t] == -(10/2) * x1[t] + (10/2) * (x2[t] - x1[t]),
        x2''[t] == -(10/2) * (x2[t] - x1[t]) - (10/2) * x2[t],
        x1[0] == 1,
        x2[0] == 0,
        x1'[0] == 0,
        x2'[0] == 0}, {x1[t], x2[t]}, t]

{{x1[t] -> 1/2 (Cos[Sqrt[5] t] + Cos[Sqrt[15] t]), 
  x2[t] -> 1/2 (Cos[Sqrt[5] t] - Cos[Sqrt[15] t])}}