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