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

Nspire で検算する。2 行目が真値、3 行目が rk23() による近似値。

TI-Nspire & Lua / 古典的ルンゲクッタ法 4 / クロージャを使って初期値を更新しながら 1 ステップずつの計算を繰り返すが、1 ステップをさらに内部分割する

function rk(func, t0, x0, step, numOfDiv) -- numOfDiv は 1 ステップの内部分割数 
   local numOfDiv = numOfDiv or 1
   local step     = step / numOfDiv
   local stepDiv2 = step / 2
   local stepDiv6 = step / 6
   local f1, f2, f3, f4

   local t0, x0    = t0, x0
   local halfPoint = t0 + stepDiv2
   local nextPoint = t0 + step
   return function()
             for i = 1, numOfDiv do
                f1 = func(t0        , x0                )
                f2 = func(halfPoint , x0 + stepDiv2 * f1)
                f3 = func(halfPoint , x0 + stepDiv2 * f2)
                f4 = func(nextPoint , x0 + step     * f3)

                t0        = nextPoint
                halfPoint = t0 + stepDiv2
                nextPoint = t0 + step
                x0        = x0 + stepDiv6 * (f1 + 2 * (f2 + f3) + f4)                
             end
             return t0, x0
          end
end

--- 確かめる。
do
local function sute(t, x)
   return x^2 - t^2 - 2 * t + 2
end

local init, from, to, step = 0, 0, 3, 1/8

local rkA = rk(sute, from, init, step, 1  ) -- 内部分割数を変えて比べてみる。t = 3 のときの真値は 3.75 である。
local rkB = rk(sute, from, init, step, 2^2)
local rkC = rk(sute, from, init, step, 2^6)

for i = 1, to/step do
   print(rkA())
   print(rkB())
   print(rkC())
end
end

TI-Nspire & Lua / 古典的ルンゲクッタ法 3 / クロージャを使って初期値を更新しながら 1 ステップずつの計算を繰り返す

function rk(func, t0, x0, step)
   local t0, x0 = t0, x0
   
   local stepDiv2  = step / 2
   local stepDiv6  = step / 6
   local halfPoint = t0 + stepDiv2
   
   local f1, f2, f3, f4

   return function()
              f1 = func(t0        , x0                )
              f2 = func(halfPoint , x0 + stepDiv2 * f1)
              f3 = func(halfPoint , x0 + stepDiv2 * f2)
              f4 = func(t0 + step , x0 + step     * f3)

              t0        = t0 + step
              x0        = x0 + stepDiv6 * (f1 + 2 * (f2 + f3) + f4)
              halfPoint = t0 + stepDiv2

              return t0, x0
           end
end

-- 確かめる。
function sute(t, x)
   return x^2 - t^2 - 2 * t + 2
end

rk1 = rk(sute, 0, 0, 1/8) -- 初期値のセットされた函数を作って、
for i = 1, 8 do
   print(rk1()) -- その初期値を更新しながら計算を繰り返す。
end