参考: パソコンで見る天体の動き, 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(listA[i] - listB[i]) end return math.max(unpack(sute)) end local dim = #funcs local f = {} local approx = {} for j = 1, dim do f[j] = funcs[j](unpack(preApprox)) end for j = 1, dim do approx[j] = w[j] + h*h * 19 * f[j]/240 -- 次の近似値 end return approx, maxOfErr(preApprox, approx) -- {次の近似値}, (前の近似値と次の近似値との最大誤差) end -- 確かめる do -- この聯立二階微分方程式で確かめる。 local function xDotDot(x, y) return -2.95912208e-4 * x / (math.sqrt(x^2 + y^2)^3) end local function yDotDot(x, y) return -2.95912208e-4 * y / (math.sqrt(x^2 + y^2)^3) end local approx2, err = approx2({xDotDot,yDotDot},{0.7231433,1.0998377}, {0.7201404,1.0952753}, 20) print("{"..table.concat(approx2," , ").."}") print(err) end