(34)二重振子 6(Nspire による二重振子シミュレーター)


local time = 0 -- シム内経過時間〔sec〕
local l1 = 1.1 --〔m〕
local l2 = 1.5 --〔m〕
local m1 = 1.7 --〔kg〕
local m2 = 2.1 --〔kg〕
local step = 0.02 --〔sec〕
local dens = 2 -- 重りの密度〔kg/m^3〕。シミュレーションとは無関係。
local vol1, vol2 = m1/dens, m2/dens -- 重りの体積〔m^3〕。シミュレーションとは無関係。
local radM1, radM2 = ((3 * vol1)/(4 * math.pi))^(1/3), ((3 * vol2)/(4 * math.pi))^(1/3) -- 重りが球体の場合の半径〔m〕。シミュレーションとは無関係。
local diaM1, diaM2 = radM1 * 2, radM2 * 2 -- 重りが球体の場合の直径〔m〕。シミュレーションとは無関係。

function on.resize()
   w, h = platform.window:width(), platform.window:height()
   x0, y0 = w/2, h/2 -- 原点を設定
   L12 = 0.95 * h/2 -- l1 + l2 の、画面上での画素数
   L1, L2 = L12 * l1/(l1 + l2), L12 * l2/(l1 + l2) -- l1 、l2 の、画面上での画素数
   unitL = L12/(l1 + l2) -- 1 メートルの、画面上での画素数
   RadM1, RadM2 = radM1 * unitL, radM2 * unitL -- 球体重りの、画面上での半径の画素数
   DiaM1, DiaM2 = diaM1 * unitL, diaM2 * unitL -- 球体重りの、画面上での直径の画素数
   fontSize = math.floor(h/17)
end

var.store("g", 9.80665)
var.store("l1", l1)
var.store("m1", m1)
var.store("l2", l2)
var.store("m2", m2)
var.store("tstart", 0)
var.store("tstop", step)
var.store("ini_theta1", 3.2) -- 初期角度 1
var.store("ini_v1", 0) -- 初速 1
var.store("ini_theta2", 3) -- 初期角度 2
var.store("ini_v2", 0)-- 初速 2
var.store("step", step)
var.store("tol", 0.0001)

function on.paint(gc)
   gc:setPen("medium")
   gc:setFont("sansserif", "r", fontSize)

   -- rk23() を計算
   mat = math.eval("calculus\\w_pend(g, l1, m1, l2, m2, tstart, tstop, ini_theta1, ini_v1, ini_theta2, ini_v2, step, tol)")
   
   -- m1 の座標を計算し、l1、m1 を描画
   x1, y1 = x0 + L1 * math.sin(mat[2][1]), y0 + L1 * math.cos(mat[2][1])
   gc:setColorRGB(0x7D7D7D)
   gc:drawLine(x0, y0, x1, y1)
   gc:drawString("l1 = "..l1.." m", fontSize, 4.2 * fontSize)
   gc:setColorRGB(0x43CD80)
   gc:fillArc(x1 - RadM1, y1 - RadM1, DiaM1, DiaM1, math.deg(mat[2][1]) + 110, 320)
   gc:drawString("m1 = "..m1.." kg", fontSize, 6.6 * fontSize)
            
   -- 支点を描画
   gc:setColorRGB(0x8FBC8F)
   gc:fillArc(x0 - 4, y0 - 4, 8, 8, 0, 360)

   -- m2 の座標を計算し、l2、m2 を描画
   x2, y2 = x1 + L2 * math.sin(mat[4][1]), y1 + L2 * math.cos(mat[4][1])
   gc:setColorRGB(0x2F4F4F)
   gc:drawLine(x1, y1, x2, y2)
   gc:drawString("l2 = "..l2.." m", fontSize, 5.4 * fontSize)
   gc:setColorRGB(0x388E8E)
   gc:fillArc(x2 - RadM2, y2 - RadM2, DiaM2, DiaM2, math.deg(mat[4][1]) + 110, 320)
   gc:drawString("m2 = "..m2.." kg", fontSize, 7.8 * fontSize)

   -- シム内の経過時間を表示
   gc:setColorRGB(0x8E8E38)
   gc:drawString("time = "..math.floor(time).." s", fontSize, fontSize)
end

function on.enterKey()
   timer.start(step)
end
function on.escapeKey()
   timer.stop()
end
function on.timer()
   platform.window:invalidate()
   
   --「次の値」を初期値に入れ換える
   var.store("ini_theta1", mat[2][2])
   var.store("ini_v1"    , mat[3][2])
   var.store("ini_theta2", mat[4][2])
   var.store("ini_v2"    , mat[5][2])
   time = time + step
end

参考:

工業力学 (機械工学基礎講座)

工業力学 (機械工学基礎講座)

(pp.229 - 231)


Excelで学ぶ微分方程式

Excelで学ぶ微分方程式

(pp.171 - 176)