参考: パソコンで見る天体の動き, pp.131-142 local function list2oneColMat(list) local oneColMat = {} for i = 1, #list do oneColMat[i] = {} oneColMat[i][1] = list[i] end return oneColMat end local function neville(...) local listOfLists = {.…

参考: パソコンで見る天体の動き, pp.131-142 h = {1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128} a = {{0.53627861}, {0.53951755}, {0.54011910}, {0.54025730}, {0.54029111}, {0.54029951}, {0.54030161}} for i = 2, #h do for j = 2, i do --a[i][j] = a[i…

参考: パソコンで見る天体の動き, pp.131-142既知のポイントの数列は前回と同じとする。 h = {4, 3, 2, 1, 1/2} a = {{277}, {-472}, {-167}, {-8}, {277/64}} for i = 2, #h do for j = 2, i do a[i][j] = a[i][j-1] + (a[i][j-1] - a[i-1][j-1]) / ((h[i-j…

参考: パソコンで見る天体の動き, pp.131-142既知のいくつかのポイントを通過する未知の曲線が偶数次多項式であると假定してその既知のポイント群の外側に存在する未知のポイントを推定する。聯立方程式を解いて求める。 この数列を既知のポイントとする。 …

参考: パソコンで見る天体の動き, pp.112-119 function cowell7(funcs, t0, inits, h) local unpack = unpack or table.unpack local t0 = t0 local inits = inits local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = m…

参考: パソコンで見る天体の動き, pp.112-119 function cowell4(funcs, t0, inits, h) local unpack = unpack or table.unpack local t0 = t0 local inits = inits local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = m…

参考: パソコンで見る天体の動き, pp.112-119 function approx2(funcs, w, preApprox, h) -- ({微分方程式}, {変化しない部分}, {前の近似値}, 刻み幅) local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = math.abs(lis…

参考: パソコンで見る天体の動き, p.112-119 function approx(funcs, w, temp, h) -- ({微分方程式}, {変化しない部分}, {假定値}, 刻み幅) local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = math.abs(listA[i] - lis…

参考: パソコンで見る天体の動き, pp.112-119 function staticVals(funcs, inits, h) --local function emptyMat(numOfRows) local sute = {} for r = 1, numOfRows do sute[r] = {} end return sute end local dim = #funcs local w = {} local f = {{},{},…

xz 面の解の軌跡を描く。y の値は圓の大きさで表現する。 function fehlberg(funcs, t0, inits, h, tol) local unpack = unpack or table.unpack local Ln = math.log local Floor = math.floor local Abs = math.abs local Max = math.max local Pow = math…

function fehlberg(funcs, t0, inits, h, tol) local unpack = unpack or table.unpack local Ln = math.log local Floor = math.floor local Abs = math.abs local Max = math.max local dim = #funcs local t0 = t0 local inits = inits local function fl…

参考: パソコンで見る天体の動き, pp.103-109 function hNew(h, err, tol) if err > tol then return 0.9 * h * (tol * h/(h * err))^(1/4) else return h end end -- 確かめる print(hNew(1, 0.0017737, 0.0001))

参考: パソコンで見る天体の動き, p.103-109 function floorB(num) if (num > 0) and (num < 1) then return 2^math.floor(math.log(num)/math.log(2)) else return math.floor(num) end end -- 確かめる print(floorB(0.1), floorB(0.2), floorB(0.8), floo…

参考: パソコンで見る天体の動き, p.103-109 local unpack = unpack or table.unpack function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = math.abs(listA[i] - listB[i]) end return math.max(unpack(sute)) end -- 確かめる…

参考: パソコンで見る天体の動き, p.103-109 function fehlberg(funcs, t0, inits, h) local unpack = unpack or table.unpack local dim = #funcs local f1 , f2 , f3 , f4 , f5 ,f6 = {}, {}, {}, {}, {}, {} local tmp1, tmp2, tmp3, tmp4, tmp5 = {}, {}…

参考: パソコンで見る天体の動き, pp.109-112特殊な方程式に対する Nyström 法を使って二階のまま解く。ここでいう「特殊な方程式」とは、一次微分も独立変数も含まない二階の方程式のことである。 function nystroem(funcs, t0, inits, initsDot, h, numOfD…

参考: Shanks, E.B., (1965), Higher Order Approximations of Runge-Kutta Type, NASA TN D-2920 (https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19650022581.pdf) function rkShanks(funcs, t0, inits, h, numOfDiv) local unpack = unpack or t…

参考: パソコンで見る天体の動き, 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 / …

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 loca…

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) …

参考: パソコンで見る天体の動き, 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 fo…

参考: パソコンで見る天体の動き, p.22 function rk(func, from, initVal, step) local f1 = func(from , initVal ) local f2 = func(from + step / 2, initVal + step * f1 / 2) local f3 = func(from + step / 2, initVal + step * f2 / 2) local f4 = fun…