（26）振子の運動 7（空気抵抗を加味する。初速をつける。）

下の連立方程式を解く。速度に比例する空気抵抗の項を -k*vel として加える。

棒の長さは 1〔m〕。
係数 k は分からないのでひとまず 0.2 とした。
初期角度は π（鉄棒に倒立しているイメージ）。
初速は 6〔rad/s〕。

横軸を角度 θ、縦軸を速度 vel（= θ'） として Nspire の Graphs アプリケーションで近似解を描いた。
よく見ると、重りは支点の真下（θ = 2π、4π、6π のところ）に来る少し手前で最速に達している。

上と同じ条件でアニメ化してみる。

local w = platform.window:width() -- 画面幅を取得
local h = platform.window:height() -- 画面高さを取得
local center = {x = w/2, y = h/2}  -- 画面中央を支点座標とする
local length = 0.95 * h/2   -- 画面上の棒の長さ
local radius = h/50 -- 画面上の重りの半径 
local diameter = radius * 2 -- 画面上の重りの直径

local g = 9.80665 --重力加速度〔m/s/s〕
local L = 1 -- 棒の長さ〔m〕
local tMin = 0 -- スタート時間 〔s〕
local tMax = 0.02 -- ストップ時間〔s〕
local angleInit = math.pi -- 角度初期値〔rad〕
   local aI = angleInit -- リセット用に角度初期値を保存しておく
local velInit = 6 -- 角速度初期値〔rad/s〕
   local vI = velInit -- リセット用に角速度初期値を保存しておく
local k = 0.2 -- 空気抵抗係数
local step = tMax -- 時間ステップ〔s〕 
local tol = 1e-4 -- 許容差
local elapsedTime = 0 -- シム内経過時間初期値

function on.paint(gc)
   --「今の角度、角速度」を初期値としてルンゲクッタ法（rk23()）で「次の角度、角速度」を計算する
   -- rk23() からは下の行列で計算結果が返ってくる。
   -- [1][1]今の時間     [1][2]次の時間
   -- [2][1]今の角度     [2][2]次の角度
   -- [3][1]今の角速度   [3][2]次の角速度
   local mat = math.eval("rk23({vel,((−"..g..")/("..L.."))*sin(θ)-"..k.."*vel},t,{θ,vel},{"..tMin..","..tMax.."},{"..angleInit..","..velInit.."},"..step..","..tol..")")
   
   -- 支点から重りの位置まで直線を引く。
   gc:drawLine(center.x, center.y, center.x + length * math.sin(mat[2][1]), center.y + length * math.cos(mat[2][1]))
   
   -- 重りの円を描く。
   gc:fillArc(center.x + length * math.sin(mat[2][1]) - radius, center.y + length * math.cos(mat[2][1]) - radius, diameter, diameter, 0, 360)   
   
   -- シム内の経過時間を表示する。
   gc:drawString("time: "..math.floor(elapsedTime).." s", 5, 0)
   elapsedTime = elapsedTime + step
         
   -- 次の角度、角速度を初期値に入れ換える。
   angleInit = mat[2][2]
   velInit = mat[3][2]
end

function on.escapeKey() -- escape キーでリセット
   angleInit = aI
   velInit = vI
   elapsedTime = 0
   timer.stop()
end
function on.enterKey() -- enter キーでタイマーをスタート
   timer.start(step)
end
function on.timer()
   platform.window:invalidate()
end

参考：

Excelで学ぶ微分方程式

（セクション 9.2.2、page 175）