(28)振子の運動 9(pendulum waves)

↓ これを Nspire + Lua + ルンゲクッタ法で再現してみる。


Pendulum Waves - YouTube


Nspire + Lua 版:


function on.resize()
   w = platform.window:width()
   h = platform.window:height()
   center = {x = w/2, y = w/80}  -- 振子の支点を設定
   radius = h/18 -- 画面上の重りの半径 
   diameter = radius * 2 -- 画面上の重りの直径

   listL = {0.864, 0.765, 0.681, 0.612, 0.553, 0.501, 0.457, 0.418, 0.384, 0.354, 0.327, 0.303, 0.282, 0.263, 0.2455} --「棒の長さ」リスト
   
   aI = 0.5 -- 角度初期値〔rad〕
   vI = 0 -- 角速度初期値〔rad/s〕   
   angleInit = math.eval("mat@>list(constructMat("..aI..",i,j,1,"..#listL.."))")
   velInit = math.eval("mat@>list(constructMat("..vI..",i,j,1,"..#listL.."))")
   
   pers =  math.eval("seq((1/n), n, 1,"..#listL..", 0.15)") -- 奥の重りほど描画上の半径を小さくして多少奥行き感を出すための係数
   --print(angleInit[1],angleInit[2])
end

local g = 9.80665 --重力加速度〔m/s/s〕
local tMin = 0 -- スタート時間 〔s〕
local tMax = 0.02 -- ストップ時間〔s〕
local k = 0.01 -- 空気抵抗係数
local step = tMax -- 時間ステップ〔s〕。1 回しか(「次の角度、角速度」しか)計算しない。
local tol = tostring(1e-3):gsub("e", string.uchar(0xF000)) -- 許容差
local elapsedTime = 0 -- シム内経過時間初期値


function on.paint(gc)
   gc:setColorRGB(0x00008B) -- 背景色を設定
   gc:fillRect(0, 0, w, h)

   mat = {}
   now = {}
   nextAngle = {}
   nextVel = {}
   for i = 15, 1, -1 do -- 奥の重り(棒の短い重り)から描画する
      gc:setColorRGB(200-(i-1)*7, 200-(i-1)*7, 255-(i-1)*2) -- 奥の重りほど色をくすませる
      -- 棒が長くなるほど弧速度(?)が速いので -(空気抵抗係数) × (棒の長さ) × (角速度) の項を加えた
      mat[i] = math.eval("rk23({vel,((−"..g..")/("..listL[i].."))*sin(θ)-"..k*listL[i].."*vel},t,{θ,vel},{"..tMin..","..tMax.."},{"..angleInit[i]..","..velInit[i].."},"..step..","..tol..")")
      now[i] = {x = center.x + h * listL[i] * math.sin(mat[i][2][1]), y = center.y + h * listL[i] * math.cos(mat[i][2][1])}
      angleInit[i] = tostring(mat[i][2][2]):gsub("e", string.uchar(0xF000))
      velInit[i] = tostring(mat[i][3][2]):gsub("e", string.uchar(0xF000))
      gc:drawLine(center.x, center.y, now[i].x, now[i].y)
      gc:fillArc(now[i].x - (radius * pers[i]), now[i].y - (radius * pers[i]), (diameter * pers[i]), (diameter * pers[i]), 0, 360) -- 奥の重りほど小さく描画する
   end
   
   -- シム内の経過時間を表示する。
   gc:setColorRGB(255, 255, 255)
   gc:setFont("sansserif", "r", h/15)
   gc:drawString("time: "..math.floor(elapsedTime).." s", w/40, h/40, "top")
   
end

---------------------------------------------------------------------------
function on.escapeKey() -- escape キーでリセット
   angleInit = math.eval("mat@>list(constructMat("..aI..",i,j,1,"..#listL.."))")
   velInit = math.eval("mat@>list(constructMat("..vI..",i,j,1,"..#listL.."))")
   elapsedTime = 0
   timer.stop()
end
function on.enterKey() -- enter キーでタイマーをスタート
   timer.start(step)
end
function on.timer()
   elapsedTime = elapsedTime + step
   platform.window:invalidate()
end


wolfram にも専用のドキュメントが用意されている。