TI-Nspire & Lua / 古典的ルンゲクッタ法 5 of 5 / 聯立方程式に適用する

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

function rkClassic(funcs, t0, inits, h, numOfDiv)
   local unpack = unpack or table.unpack  

   local dim      = #funcs
   local numOfDiv = numOfDiv or 1

   local h     = h / numOfDiv
   local hDiv2 = h / 2
   local hDiv6 = h / 6

   local f1, f2, f3, f4         = {}, {}, {}, {}
   local temps1, temps2, temps3 = {}, {}, {}

   local t0        = t0
   local inits     = inits
   local halfPoint = t0 + hDiv2
   local nextPoint = t0 + h
   return function() 
             for i = 1, numOfDiv do   
                for j = 1, dim do f1[j] = funcs[j](t0       , unpack(inits))  ; temps1[j] = inits[j] + hDiv2 * f1[j] end
                for j = 1, dim do f2[j] = funcs[j](halfPoint, unpack(temps1)) ; temps2[j] = inits[j] + hDiv2 * f2[j] end
                for j = 1, dim do f3[j] = funcs[j](halfPoint, unpack(temps2)) ; temps3[j] = inits[j] + h     * f3[j] end
                for j = 1, dim do f4[j] = funcs[j](nextPoint, unpack(temps3))                                        end
                
                for j = 1, dim do inits[j] = inits[j] + hDiv6 * (f1[j] + 2 * (f2[j] + f3[j]) + f4[j])                end 

                t0        = nextPoint
                halfPoint = t0 + hDiv2
                nextPoint = t0 + h 
             end
             return t0, unpack(inits)
          end
end


--- 確かめる。
do

-- この聯立微分方程式を解く。
local function xDot(t, x, y) return y     end -- x'(t) = y(t)     
local function yDot(t, x, y) return t - x end -- y'(t) = t - x(t)

-- 1 ステップだけ計算する函数を作って、
local rk = rkClassic({xDot, yDot}, 0, {0, 0}, 0.25, 1) -- (函数リスト, t0, 初期値リスト, 刻み幅 [, 内部分割数])

-- その函数を必要な回数だけ繰り返す。
for i = 1, 7.5/0.25 do
   print(rk())
end

end

f:id:ti-nspire:20170505073530p:plain:h300


Nspire で検算する。2 行目が真値、3 行目が rk23() による近似値。
f:id:ti-nspire:20170505073713p:plain