TI-Nspire & Lua / 古典的ルンゲクッタ法 2 / 何ステップかを一度にまとめて計算する

参考: パソコンで見る天体の動き, p.25

function rk(func, period, init, step)
   local step     = step
   local from     = period[1]
   local to       = period[2]
   local initList = {init}
   local timeList = {from}
   local stepDiv6 = step / 6
   local stepDiv2 = step / 2

   for t = from, (to - step) + (step / 2), step do
      local init      = initList[#initList]
      local halfPoint = t + stepDiv2
      local nextPoint = t + step

      local f1 = func(t         , init                )
      local f2 = func(halfPoint , init + stepDiv2 * f1)
      local f3 = func(halfPoint , init + stepDiv2 * f2)
      local f4 = func(nextPoint , init + step     * f3)

      initList[#initList+1] = init + stepDiv6 * (f1 + 2 * (f2 + f3) + f4)
      timeList[#timeList+1] = nextPoint
   end
   var.store("rk", {})
   var.store("rk", {timeList, initList})
end


-- 確かめる。
function sute(t, x)
   return x^2 - t^2 - 2 * t + 2
end
rk(sute, {0, 1}, 0, 1/8)

Nspire で検算。2 行目が厳密解、3 行目が rk23() による近似解、4 行目が今回の古典ルンゲクッタ法による近似解。