TI-Nspire & Lua / 補外法 3 / 何かを偶数次多項式で補外する / Neville の算法で補外値だけを求める / 刻み幅を順次 1/2 にした場合

参考: パソコンで見る天体の動き, 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][j-1] + (a[i][j-1] - a[i-1][j-1]) / ((h[i-j+1] / h[i])^2 - 1)
   a[i][j] = a[i][j-1] + (a[i][j-1] - a[i-1][j-1]) / (4^(j-1) - 1) -- 刻み幅を順次 1/2 にした場合は公式が少し単純になる。
end
end

-- 確かめる。
for i = 1, #a do
   print("{"..table.concat(a[i], ", ").."}")
end

右下の 0.54030231009926 が補外値。

―――――――――――――――――――――――――――――――――
Nspire で聯立方程式を解いて検算する。

x:={((1)/(2)),((1)/(4)),((1)/(8)),((1)/(16)),((1)/(32)),((1)/(64)),((1)/(128))}
y:={0.53627861,0.53951755,0.5401191,0.5402573,0.54029111,0.54029951,0.54030161}
linSolve(seqn(y[n]=polyEval({a6,a5,a4,a3,a2,a1,a0},x[n]^(2)),7),a6,a5,a4,a3,a2,a1,a0)

右端の 0.540302310099 が補外値。