Adaで関数の近似(ネヴィル補間)

をネヴィル補間で近似する

与えられた５個の関数値 f(x₀), f(x₁), f(x₂), f(x₃), f(x₄) を通る４次式は、次のようにして求めることができる.
(1)

(2)

(3)

(4)

(5)

これらの式を使って, 与えられた点以外の点の値を求める.

with Text_IO, Ada.Long_Float_Text_IO;
use  Text_IO, Ada.Long_Float_Text_IO;

procedure Ada0702 is
    N : Constant Integer := 7;
    type Float_Vector is array (0..N-1) of Long_Float;
    x : Float_Vector;
    y : Float_Vector;
    d, d1, d2 : 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;
    
    -- Neville (ネヴィル) 補間
    function neville(d:Long_Float; x:Float_Vector; y:Float_Vector) return Long_Float is
        w : array (0..N-1, 0..N-1) of Long_Float;
    begin
        for i in y'Range loop
            w(0, i) := y(i);
        end loop;

        for j in 1 .. x'Last loop
            for i in 0 .. x'Last-j loop
                w(j,i) := w(j-1,i+1) + (w(j-1,i+1) - w(j-1,i)) * (d - x(i+j)) / (x(i+j) - x(i));
            end loop;
        end loop;
        return w(N-1,0);
    end;
begin
    -- 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for i in 0..N-1 loop
        d    := Long_Float(i - 1) * 1.5 - 4.5;
        x(i) := d;
        y(i) := f(d);
    end loop;
    
    -- 0.1刻みで 与えられていない値を補間 
    for i in 0..90 loop
        d  := Long_Float(i) * 0.1 - 4.5;
        d1 := f(d);
        d2 := neville(d, x, y);

        -- 元の関数と比較
        Put(d,       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 Ada0702;

Compiling the source code....
$gnatmake ada0702.adb 2>&1
gcc -c ada0702.adb
gnatbind -x ada0702.ali
gnatlink ada0702.ali

Executing the program....
$ada0702 
-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
以下略

参考文献