TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 2 / 第 1 近似値、第 2 近似値、両者の最大誤差を求める、それに意味のないとき

参考: パソコンで見る天体の動き, 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] - listB[i])
            end
            return math.max(unpack(sute))
         end

   local dim     = #funcs
   local approx1 = {}
   local approx2 = {}
   local f       = {}

   for j = 1, dim do
      approx1[j] = w[j] + h*h * 19 * temp[j]/240 -- 第 1 近似値
   end
   for j = 1, dim do
      f[j] = funcs[j](unpack(approx1)) 
   end
   for j = 1, dim do
         approx2[j] = w[j] + h*h * 19 * f[j]/240 -- 第 2 近似値
   end
   
   return approx1, approx2, maxOfErr(approx1, approx2) -- {第 1 近似値}, {第 2 近似値}, 最大誤差  
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 approx1, approx2, err = approx({xDotDot,yDotDot},{0.7231433,1.0998377},{-0.9547563e-4,-1.8131274e-4},20)
print("{"..table.concat(approx1," , ").."}")
print("{"..table.concat(approx2," , ").."}")
print(err)

end