Adaで関数の近似(ニュートン補間)

をニュートン補間で近似する

与えられた４個の関数値 f(x₀), f(x₁), f(x₂), f(x₃) を通る３次式を求める場合, まず次のような表を作る.

このとき,

と定義する.
これを x₀とx₁ の第１差分商, または１階差分商という.
同様に

と定義し, これを x₀, x₁, x₂ の第２差分商, または２階差分商という.
第ｎ差分商を

で表すと, 与えられた 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 まで, ７点だけ値をセット
    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;

    -- ｎ階差分商
    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
以下略

ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ立ったら歩きなはれ

参考文献