逆数補間の計算方法について

質問

こんにちは。前にも書かせてもらいましたが、どうしても計算ができないので、もう一度質問させてもらいました。
以下のような、洋書を読んで、最後にあるP(y)を出したいのですが、計算方法がわかりません。
----------------------------------------------------------------
[Inverse Interpolation]
A process called inverse interpolation is often used to approximate an inverse function. Suppose that values {Yi}=f({Xi}) have been computed at X0,X1,...,Xn.
Using table
Y ; Y0 Y1 Y2 ......Yn
X ; X0 X1 X2 ......Xn
we form the interpolation polynomial

p(y)=Σ(i=1→n)CiΠ(j=0→i-1){Y－Yj}
The orijinal relationship, y=f(x), has an inverse, under certain conditions. This inverse is being approximated by x=p(y). Procedures Coef and Eval can be used to carry out the inverse interpolation by reversing the arguments x and y in the calling sequence for Coef.
Inverse interpolation can be used to find where a given functuin f has a root or zero. This means inverting the equation f(x)=0. We propose to do this by creating a table of values (f(Xi),Xi) and interpolating with a polynomial,p. Thus, p(Yi)=Xi. The points Xi should be chosen near the unknown root,r. The approximate root is then given by r ～p(0). For a concrete case, let the table of known values be

Y;-0.5789200,-0.3626370,-0.1849160,-0.0340642,0.0969858

X; 1.0 , 2.0 , 3.0 , 4.0 , 5.0
The nodes in this problem are the points in the row of the table headed y, and the function values being interpolated are in the x row. The resulting polynomial is

p(Y)=0.25Y^4＋1.2Y^3＋3.69Y^2＋7.39Y＋4.247470086

and p(0)=4.247470086. Only the last coefficient is shown with all the digits carried in the calculation, for it is the only one needed for the problem at hand.
----------------------------------------------------------------
＜補足＞CoefとEvalについて

「 procedure; Coef(n,{Xi},{Yi},{Ai})
real array; {Xi}0:n, {Yi}0:n, {Ai}0:n
integer; i,j,n
for i=0 to n do
{Ai}←{Yi}
end for
for j=1 to n do
for i=n to j step -1 do
Ai←({Ai}－{Ai-1})/({Xi}－{Xi-j})
end for
end for
end procedure Coef 」

「 real function; Eval(n,{Xi},{Yi},{Ai})
real array; {Xi}0:n, {Ai}0:n
integer; i,n
real;t,temp
temp←An
for i=n-1 to 0 step -1 do
temp←(temp)(t-{Xi})＋{Ai}
end for
Eval←temp
end function Eval」
-------------------------------------------------------------
ＸとＹを扱い方がよくわかっていないので、計算できないのかなあと思います。分かる方、アドバイスお願いします（泣）

回答No.4

noname#101087 · ベストアンサー

切れぎれですみません。

>Y;-0.5789200,-0.3626370,-0.1849160,-0.0340642,0.0969858
>X; 1.0 , 2.0 , 3.0 , 4.0 , 5.0
> .........
>p(Y)=0.25Y^4＋1.2Y^3＋3.69Y^2＋7.39Y＋4.247470086
>----------------------------------------------------------------
>＜補足＞CoefとEvalについて
>..........
>-------------------------------------------------------------
>ＸとＹを扱い方がよくわかっていないので、計算できないのかなあと思います。....

Coef と Eval は良さそうにみえます。

ただ、Coef で求めているのは、
　P4(x) = A0 + A1(x - x0) + A2(x - x0)(x - x1) + A3(x - x0)(x - x1)(x - x2) + A4(x - x0)(x - x1)(x - x2)(x - x3)
の係数{A0, A1, .... , A4}ですね。(Eval を見ても明らかです)

つまり、[Inverse Interpolation]の例である
　p(Y)=0.25Y^4＋1.2Y^3＋3.69Y^2＋7.39Y＋4.247470086
の係数が Coef の結果と違うのは当然なのです。
　(1) [Inverse Interpolation]に例示された対応表 {Xi,Yi} の X,Y を入れ替え
　(2) 上記の P4(x) を算出
　(3) P4(x) をふつうの多項式に展開
すれば、引用の p(Y) になるでしょう。

そこまでやらずとも、Coef で求めた係数{A0, A1, .... , A4} を使い Eval で計算してみれば、対応表 {Xi,Yi} と合っているか否か
はわかるはずですが....。

回答No.3

noname#101087

>procedure; Coef(n, {Xi}, {Yi}, {Ai})
>......
>for i=0 to n do
>{Ai}←{Yi}
>end for
>for j=1 to n do
>for i=n to j step -1 do
>Ai←({Ai}－{Ai-1})/({Xi}－{Xi-j})
>end for
>......
-------------------
参考までに Coef における{Ai}(差分商表)の書き換えをスケッチしてみます。
「多分」OK だと思います。

(1) 係数は、小文字(ai)が計算途中、大文字(Ai)が計算結果、と区別。
(2) Xn-Xm=dnm と略記。

j =0　　j =1　　　　　　　　　　　　　j =2
---　　---　　　　　　　　　　　　　---
Y0 　Y0=A0
Y1 　(Y1-Y0)/d10=A1
Y2 　(Y2-Y1)/d21=a2　　　　　　　(a2-A1)/d20=A2
.......
Yn 　(Yn-Y<n-1>)/dn<n-1>=an

回答No.2

noname#101087

Newton 補間の書き忘れです。

下記ページなど、参考になりませんか?
------------------------------------
　http://www.activebasic.com/aba/wiki.cgi?Newton%28%a5%cb%a5%e5%a1%bc%a5%c8%a5%f3%29%ca%e4%b4%d6
>Newton(ニュートン)補間

回答No.1

noname#101087

>最後にあるP(y)を出したいのですが、計算方法がわかりません。

引っかかりポイントか P(Y) の作り方らしいので、まとめてみましょう。
(QNo.2980931 & QNo.2980811 の繰り返しですが.... )

(0) 補間表。
　　Y0=-0.5789200, Y1=-0.3626370, Y2=-0.1849160, Y3=-0.0340642, Y4=0.0969858
　　X0=1.0, X1=2.0, X2=3.0, X3=4.0, X4=5.0

(1) 補間多項式の作成。
　　X=P4(Y)=0.2504105*Y^4+1.2156411*Y^3+3.6940309*Y^2+7.3892454*Y+4.2474700

　当方は、EXCEL シートで作りました。
　Lagrange 補間公式をそのまま使ってます。(記号 X, Y は入れ替わってます)
　一例。Y0 にてX0 になり、ほかの{X1, X2, X3, X4}では零になる多項式は
　　Q0(Y)=A0*(X-X1)*(X-X2)*(X-X3)*(X-X4)
　です。A0 は、Q0(Y0)=X0 となるように算定します。
　同様に、Q1(Y), Q2(Y), Q3(Y), Q4(Y) を求め、すべてを合計したのが上記の P(Y) です。
　EXCEL シートだと、Q0(Y)用の式をコピーして{Q1(Y), Q2(Y), Q3(Y), Q4(Y)} を楽に作れるのです。
　(EXCEL シートは、一覧性が良くミスを見つけ易い)

　シーケンシャルなプログラムの場合は、Newton 補間公式のほうがコンパクトになるでしよう。(結果は同じになります)