ONLY DO WHAT ONLY YOU CAN DO

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

Adaで関数の近似(ラグランジュ補間)



ラグランジュ補間で近似する

n+1個の点 (x0, y0), (x1, y1) … (xn, yn)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.

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

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

procedure Ada0701 is
    N : Constant Integer := 7;
    type Float_Vector is array (1..N) 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;
    
    -- Lagrange (ラグランジュ) 補間
    function lagrange(d:Long_Float; x:Float_Vector; y:Float_Vector) return Long_Float is
        sum1, prod :Long_Float;
    begin
        sum1 := Long_Float(0.0);
        --for i in x'First .. x'Last loop
        for i in x'Range loop
            prod := y(i);

            for j in x'Range loop
            --for j in x'First .. x'Last loop
                if j /= i then
                    prod := prod * (d - x(j)) / (x(i) - x(j));
                end if;
            end loop;
            sum1 := sum1 + prod;
        end loop;
        return sum1;
    end;
begin
    -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット
    for i in 1..N 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 := lagrange(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 Ada0701;
Compiling the source code....
$gnatmake ada0701.adb 2>&1
gcc -c ada0701.adb
gnatbind -x ada0701.ali
gnatlink ada0701.ali

Executing the program....
$ada0701 
-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
-1.90	 -0.96317	 -0.96317	 -0.00000
-1.80	 -0.98546	 -0.98546	 -0.00000
-1.70	 -0.99949	 -0.99949	 -0.00000
-1.60	 -1.00471	 -1.00471	 -0.00000
-1.50	 -1.00078	 -1.00078	  0.00000
-1.40	 -0.98749	 -0.98749	  0.00000
-1.30	 -0.96477	 -0.96477	  0.00000
-1.20	 -0.93274	 -0.93274	 -0.00000
-1.10	 -0.89159	 -0.89159	  0.00000
-1.00	 -0.84167	 -0.84167	  0.00000
-0.90	 -0.78342	 -0.78342	  0.00000
-0.80	 -0.71740	 -0.71740	  0.00000
-0.70	 -0.64423	 -0.64423	  0.00000
-0.60	 -0.56465	 -0.56465	  0.00000
-0.50	 -0.47943	 -0.47943	 -0.00000
-0.40	 -0.38942	 -0.38942	 -0.00000
-0.30	 -0.29552	 -0.29552	  0.00000
-0.20	 -0.19867	 -0.19867	 -0.00000
-0.10	 -0.09983	 -0.09983	  0.00000
 0.00	  0.00000	  0.00000	  0.00000
 0.10	  0.09983	  0.09983	  0.00000
 0.20	  0.19867	  0.19867	  0.00000
 0.30	  0.29552	  0.29552	 -0.00000
 0.40	  0.38942	  0.38942	 -0.00000
 0.50	  0.47943	  0.47943	  0.00000
 0.60	  0.56465	  0.56465	  0.00000
 0.70	  0.64423	  0.64423	  0.00000
 0.80	  0.71740	  0.71740	 -0.00000
 0.90	  0.78342	  0.78342	  0.00000
 1.00	  0.84167	  0.84167	  0.00000
 1.10	  0.89159	  0.89159	  0.00000
 1.20	  0.93274	  0.93274	  0.00000
 1.30	  0.96477	  0.96477	  0.00000
 1.40	  0.98749	  0.98749	  0.00000
 1.50	  1.00078	  1.00078	  0.00000
 1.60	  1.00471	  1.00471	 -0.00000
 1.70	  0.99949	  0.99949	  0.00000
 1.80	  0.98546	  0.98546	  0.00000
 1.90	  0.96317	  0.96317	 -0.00000
 2.00	  0.93333	  0.93333	 -0.00000
 2.10	  0.89684	  0.89684	  0.00000
 2.20	  0.85480	  0.85480	  0.00000
 2.30	  0.80853	  0.80853	 -0.00000
 2.40	  0.75955	  0.75955	  0.00000
 2.50	  0.70964	  0.70964	  0.00000
 2.60	  0.66078	  0.66078	 -0.00000
 2.70	  0.61524	  0.61524	  0.00000
 2.80	  0.57553	  0.57553	  0.00000
 2.90	  0.54443	  0.54443	 -0.00000
 3.00	  0.52500	  0.52500	  0.00000
 3.10	  0.52060	  0.52060	  0.00000
 3.20	  0.53487	  0.53487	 -0.00000
 3.30	  0.57178	  0.57178	  0.00000
 3.40	  0.63562	  0.63562	  0.00000
 3.50	  0.73099	  0.73099	  0.00000
 3.60	  0.86285	  0.86285	  0.00000
 3.70	  1.03650	  1.03650	  0.00000
 3.80	  1.25760	  1.25760	  0.00000
 3.90	  1.53218	  1.53218	  0.00000
 4.00	  1.86667	  1.86667	  0.00000
 4.10	  2.26785	  2.26785	 -0.00000
 4.20	  2.74294	  2.74294	  0.00000
 4.30	  3.29954	  3.29954	  0.00000
 4.40	  3.94569	  3.94569	 -0.00000
 4.50	  4.68984	  4.68984	  0.00000
参考文献