TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 重振子でバタフライ効果を見る / 古典的 Runge-Kutta 法

L2 の長さを 0.01% だけ違えた 2 組の 2 重振子を同一座標に描く。t = 6 あたりまでは同じ軌道を描くが、t = 10 を過ぎると全然違う軌道になる。

DoublePend = class()
function DoublePend:init(funcs, t0, inits, h, numOfDiv, len1, radius1, len2, radius2)
   self.len1    = len1
   self.radius1 = radius1

   self.len2    = len2
   self.radius2 = radius2
   
   self.historicData = {{},{}}
  
   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 = {{}, {}, {}}
   
   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
end
function DoublePend: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 DoublePend:print()
   print(self.t0, "{"..table.concat(self.inits, ", ").."}")
   return self
end
function DoublePend: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 DoublePend:drawPend(gc, cx, cy, unit, color)
   local len1OnScr    = unit * self.len1
   local radius1OnScr = unit * self.radius1
   local dia1OnScr    = 2 * radius1OnScr

   local len2OnScr    = unit * self.len2
   local radius2OnScr = unit * self.radius2
   local dia2OnScr    = 2 * radius2OnScr
   
   self.sx    = {}
   self.sy    = {}
   self.sx[1] = cx         + len1OnScr * self.Sin(self.inits[1])
   self.sy[1] = cy         + len1OnScr * self.Cos(self.inits[1])
   self.sx[2] = self.sx[1] + len2OnScr * self.Sin(self.inits[3])
   self.sy[2] = self.sy[1] + len2OnScr * self.Cos(self.inits[3])

   gc:setColorRGB(0x000000)
   gc:setPen("thin")
   gc:drawLine(cx, cy, self.sx[1], self.sy[1])
   gc:drawLine(self.sx[1], self.sy[1], self.sx[2], self.sy[2])

   gc:setColorRGB(color or 0x000000)
   gc:setPen("medium")
   gc:drawArc(self.sx[2] - radius2OnScr, self.sy[2] - radius2OnScr, dia2OnScr, dia2OnScr, 0, 360)
   gc:drawArc(self.sx[1] - radius1OnScr, self.sy[1] - radius1OnScr, dia1OnScr, dia1OnScr, 0, 360)
   
   return self
end
function DoublePend:drawLocus(gc, which, numOfPairs, color)
   self.storeData(self.sx[which], self.sy[which], self.historicData[which], numOfPairs)
   
   gc:setPen("thin")
   gc:setColorRGB(color or 0x000000)
   gc:drawPolyLine(self.historicData[which])
   
   return self
end


--- 確かめる。
do
require "color"
platform.window:setBackgroundColor(0xE7E9DE)
local colorsPend  = {0x5D86FF, 0x78B9FF}
local colorsLocus = {0x7D9EC0, 0x6CA6CD}

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/4.5
   FONT_SIZE = H/16
end

local Sin = math.sin  
local Cos = math.cos

-- この聯立微分方程式を解く。
local function rad1Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return vel1 end
local function vel1Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return -(9.80665*(m1*Sin(rad1)-m2*(Cos(rad1-rad2)*Sin(rad2)-Sin(rad1)))+(l1*vel1^2*Cos(rad1-rad2)+l2*vel2^2)*m2*Sin(rad1-rad2))/(l1*(m1+m2*(Sin(rad1-rad2))^2)) end
local function rad2Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return vel2 end
local function vel2Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return (9.80665*(m1+m2)*(Sin(rad1)*Cos(rad1-rad2)-Sin(rad2))+(l1*(m1+m2)*vel1^2+l2*m2*vel2^2*Cos(rad1-rad2))*Sin(rad1-rad2))/(l2*(m1+m2*(Sin(rad1-rad2))^2)) end
local function len1Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return 0 end
local function len2Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return 0 end
local function mas1Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return 0 end
local function mas2Dot(t, rad1, vel1, rad2, vel2, m1, l1, m2, l2) return 0 end
local funcs = {rad1Dot, vel1Dot, rad2Dot, vel2Dot, len1Dot, len2Dot, mas1Dot, mas2Dot}

-- この初期値で解く。
local mass1, mass2 = {}, {}
mass1[1], mass2[1] = 1, 1
mass1[2], mass2[2] = 1, 1

local len1init, len2init = {}, {}
len1init[1], len2init[1] = 1, 1
len1init[2], len2init[2] = 1, 1.0001 -- L2 の長さを 0.01% だけ違える。 

local rad1init, rad2init = {}, {}
rad1init[1], rad2init[1] = 3, 3
rad1init[2], rad2init[2] = 3, 3

local vel1init, vel2init = {}, {}
vel1init[1], vel2init[1] = 0, 0
vel1init[2], vel2init[2] = 0, 0

local radius = 0.2


-- インスタンスを 2 個作る。
local inits = {}
local pends = {}
local function reset()
         for i = 1, 2 do
            inits[i] = {rad1init[i], vel1init[i], rad2init[i], vel2init[i], mass1[i], len1init[i], mass2[i], len2init[i]}
            pends[i] = DoublePend(funcs, 0, inits[i], 1/64, _, len1init[i], radius, len2init[i], radius)
         end
      end


function on.construction()
   reset()
end
function on.enterKey()
   if     isTicking == false then timer.start(pends[1].h * pends[1].numOfDiv); isTicking = true
   elseif isTicking == true  then timer.stop()                               ; isTicking = false
   end
end
function on.escapeKey()
   reset()
end
function on.timer()
   for i = 1, #pends do
      pends[i]:updateRK()
   end
   platform.window:invalidate()
end
function on.paint(gc)
   pends[1]:drawTime(gc, 10, 0, FONT_SIZE, 0x444444)
   for i = 1, #pends do
      pends[i]:drawPend(gc, CX, CY, UNIT, colorsPend[i])
              :drawLocus(gc, 2, 50, colorsLocus[i]) -- 質点 2 のスクリーン座標を過去 50 ペアぶん描く。
   end
end

end