Skip to content

Instantly share code, notes, and snippets.

@ti-nspire
Last active December 20, 2017 01:43
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save ti-nspire/7348f03ac601b7626addc122c0290b90 to your computer and use it in GitHub Desktop.
rkClassicClass.lua
-- Lua & TI-Nspire による常微分方程式の数値解法 / 古典的ルンゲクッタ法
rkClassic = class()
function rkClassic:init(funcs, t0, inits, h, numOfDiv)
self.Unpack = unpack or table.unpack
self.funcs = funcs
self.t0 = t0
self.inits = inits
self.dim = #funcs
self.numOfDiv = numOfDiv or 1
self.h = h / self.numOfDiv
self.f = {{}, {}, {}, {}}
self.temp = {{}, {}, {}}
end
function rkClassic:updateRK()
for i = 1, self.numOfDiv do
for j = 1, self.dim do self.f[1][j] = self.funcs[j](self.t0 , self.Unpack(self.inits)) ; self.temp[1][j] = self.inits[j] + self.h/2 * self.f[1][j] end
for j = 1, self.dim do self.f[2][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[1])) ; self.temp[2][j] = self.inits[j] + self.h/2 * self.f[2][j] end
for j = 1, self.dim do self.f[3][j] = self.funcs[j](self.t0 + self.h/2, self.Unpack(self.temp[2])) ; self.temp[3][j] = self.inits[j] + self.h * self.f[3][j] end
for j = 1, self.dim do self.f[4][j] = self.funcs[j](self.t0 + self.h , self.Unpack(self.temp[3])) end
for j = 1, self.dim do self.inits[j] = self.inits[j] + self.h/6 * (self.f[1][j] + 2 * (self.f[2][j] + self.f[3][j]) + self.f[4][j]) end
self.t0 = self.t0 + self.h
end
return self
end
function rkClassic:print()
print(self.t0, "{"..table.concat(self.inits, ", ").."}")
return self
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 rk = rkClassic({xDot, yDot}, 0, {0, 0}, 0.25, _) -- (函数リスト, t0, 初期値リスト, 刻み幅 [, 内部分割数])
-- 必要な回数だけそのインスタンスを更新し、更新するたびに print して確認する。
for i = 1, 30 do
rk:updateRK():print()
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment