Adaで関数の近似(ニュートン補間)
をニュートン補間で近似する
与えられた4個の関数値 f(x0), f(x1), f(x2), f(x3) を通る3次式を求める場合, まず次のような表を作る.
このとき,
と定義する.
これを x0とx1 の第1差分商, または1階差分商という.
同様に
と定義し, これを x0, x1, x2 の第2差分商, または2階差分商という.
第n差分商を
で表すと, 与えられた n+1点を通る n次式は次のように表すことができる.
この式を使って, 与えられた点以外の点の値を求める.
with Text_IO, Ada.Long_Float_Text_IO; use Text_IO, Ada.Long_Float_Text_IO; procedure Ada0703 is N : Constant Integer := 7; type Float_Vector is array (0..N-1) of Long_Float; x : Float_Vector; y : Float_Vector; a : Float_Vector; d : array (0..N-1, 0..N-1) of Long_Float; d1, d2, d3 : Long_Float; -- 元の関数 function f(x:Long_Float) return Long_Float is begin return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2); end f; -- Newton (ニュートン) 補間 function newton(d:Long_Float; x:Float_Vector; a:Float_Vector) return Long_Float is sum, prod :Long_Float; begin sum := a(0); for i in 1 .. x'Last loop prod := 1.0; for j in x'First .. i-1 loop prod := prod * (d - x(j)); end loop; sum := sum + (a(i) * prod); end loop; return sum; end; begin -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i in 0..N-1 loop d3 := Long_Float(i - 1) * 1.5 - 4.5; x(i) := d3; y(i) := f(d3); end loop; -- 差分商の表を作る for j in 0 .. N-1 loop d(0,j) := y(j); end loop; for i in 1 .. N-1 loop for j in 0 .. N-i-1 loop d(i,j) := (d(i-1,j+1) - d(i-1,j)) / (x(j+i) - x(j)); end loop; end loop; -- n階差分商 for j in 0 .. N-1 loop a(j) := d(j,0); end loop; -- 0.1刻みで 与えられていない値を補間 for i in 0..90 loop d3 := Long_Float(i) * 0.1 - 4.5; d1 := f(d3); d2 := newton(d3, x, a); -- 元の関数と比較 Put(d3, Fore=>2, Aft=>2, Exp=>0); Put(Ascii.HT); Put(d1, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d2, Fore=>3, Aft=>5, Exp=>0); Put(Ascii.HT); Put(d1 - d2, Fore=>3, Aft=>5, Exp=>0); New_Line; end loop; end Ada0703;
Compiling the source code.... $gnatmake ada0703.adb 2>&1 gcc -c ada0703.adb gnatbind -x ada0703.ali gnatlink ada0703.ali Executing the program.... $ada0703 -4.50 -4.68984 -4.68984 0.00000 -4.40 -3.94569 -3.94569 0.00000 -4.30 -3.29954 -3.29954 0.00000 -4.20 -2.74294 -2.74294 -0.00000 -4.10 -2.26785 -2.26785 -0.00000 -4.00 -1.86667 -1.86667 0.00000 -3.90 -1.53218 -1.53218 -0.00000 -3.80 -1.25760 -1.25760 0.00000 -3.70 -1.03650 -1.03650 -0.00000 -3.60 -0.86285 -0.86285 0.00000 -3.50 -0.73099 -0.73099 0.00000 -3.40 -0.63562 -0.63562 -0.00000 -3.30 -0.57178 -0.57178 0.00000 -3.20 -0.53487 -0.53487 -0.00000 -3.10 -0.52060 -0.52060 -0.00000 -3.00 -0.52500 -0.52500 0.00000 -2.90 -0.54443 -0.54443 0.00000 -2.80 -0.57553 -0.57553 -0.00000 -2.70 -0.61524 -0.61524 0.00000 -2.60 -0.66078 -0.66078 0.00000 -2.50 -0.70964 -0.70964 0.00000 -2.40 -0.75955 -0.75955 0.00000 -2.30 -0.80853 -0.80853 0.00000 -2.20 -0.85480 -0.85480 0.00000 -2.10 -0.89684 -0.89684 0.00000 -2.00 -0.93333 -0.93333 0.00000 以下略