rkClassic.lua

rkClassic.lua

-- Lua で常微分方程式を解く / 古典的ルンゲクッタ法
-- http://ti-nspire.hatenablog.com/entry/2017/07/19/143509
function rkClassic(funcs, t0, inits, h, numOfDiv)
   local unpack = unpack or table.unpack  

   local step  = h
   local t0    = t0
   local inits = inits

   local dim      = #funcs
   local numOfDiv = numOfDiv or 1
   local h        = h / numOfDiv

   local f    = {{}, {}, {}, {}}
   local temp = {{}, {}, {}}

   return function() 
             for i = 1, numOfDiv do   
                for j = 1, dim do f[1][j] = funcs[j](t0      , unpack(inits))   ; temp[1][j] = inits[j] + h/2 * f[1][j] end
                for j = 1, dim do f[2][j] = funcs[j](t0 + h/2, unpack(temp[1])) ; temp[2][j] = inits[j] + h/2 * f[2][j] end
                for j = 1, dim do f[3][j] = funcs[j](t0 + h/2, unpack(temp[2])) ; temp[3][j] = inits[j] + h   * f[3][j] end
                for j = 1, dim do f[4][j] = funcs[j](t0 + h  , unpack(temp[3]))                                         end
                
                for j = 1, dim do inits[j] = inits[j] + h/6 * (f[1][j] + 2 * (f[2][j] + f[3][j]) + f[4][j]) end 

                t0 = t0 + h
             end
             return t0, inits -- 更新した t0, {更新した初期値} という形で返す。
          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)
local funcs = {xDot, yDot}

-- この初期値で解く。
local x0 = 0
local y0 = 0
local inits = {x0, y0}

-- この刻み幅で計算する。
local h = 1/2^2

-- 独立変数の変化幅は下のとおりとする。
local t0   = 0
local tMax = 7.5

-- 1 ステップだけ計算する函数を作って、
local a = rkClassic(funcs, t0, inits, h, _) -- (函数リスト, t0, 初期値リスト, 刻み幅 [, 内部分割数])

-- その函数を必要な回数だけ繰り返す。
for i = 1, tMax/h do
   local t0, inits = a() 
   print(t0, table.concat(inits, "\t"))
end

end

rkClassic_results.lua

0.25	0.0026041666666667	0.031087239583333
0.5	0.020590040418837	0.12241276105245
0.75	0.068382146388844	0.26829855226377
1.0	0.15855187449442	0.45967454738203
1.25	0.30103590834733	0.68464253439676
1.5	0.50251845941971	0.92921588974218
1.75	0.76601571310769	1.1781891493029
2.0	1.0906883214633	1.41608335317
2.25	1.4718935993635	1.6281083907547
2.5	1.9014741782159	1.8010825161345
2.75	2.3682651704823	1.9242518664812
3.0	2.8587883133029	1.9899590335334
3.25	3.358089934906	1.9941191227871
3.5	3.8506706460176	1.9364737039217
3.75	4.3214489560065	1.8206068648273
4.0	4.7566989045644	1.653722371938
4.25	5.1449034154	1.4461957918634
4.5	5.4774703182794	1.210929420175
4.75	5.7492665239882	0.96255012280547
5.0	5.9589371425726	0.71649996168602
5.25	6.1089897056474	0.48807614179793
5.5	6.2056382573652	0.2914799673729
5.75	6.2584170078942	0.13893393466403
6.0	6.2795875697727	0.039921852385396
6.25	6.2833766306124	0.00059923196546536
6.5	6.2850914574582	0.023410603736914
6.75	6.3001682232957	0.10693755024568
7.0	6.3432123225346	0.24598690384165
7.25	6.427090340151	0.43191362440051
7.5	6.5621301276054	0.65315828271919