参考: 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
- 作者:鈴木 肇
- 発売日: 2006/12/01
- メディア: 単行本
―――――――――――――――――――
厳密解が求まるので近似解を求める意味は計算の確認以外にはない。上と同じ条件で厳密解を求めると下のようになる。
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])}}