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
参考文献