Pascalで関数の近似(ラグランジュ補間)
をラグランジュ補間で近似する
n+1個の点 (x0, y0), (x1, y1) … (xn, yn)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.
この式を使って, 与えられた点以外の点の値を求める.
Program Pas0701(arg); {$MODE delphi} uses SysUtils, Math; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Lagrange (ラグランジュ) 補間 function lagrange(d:Double; x:array of Double; y:array of Double):Double; var sum, prod :Double; i, j :Integer; begin sum := 0; for i := Low(x) to High(x) do begin prod := y[i]; for j := Low(x) to High(x) do begin if j <> i then prod := prod * ((d - x[j]) / (x[i] - x[j])); end; sum := sum + prod; end; result := sum; end; const // データ点の数 N = 7; var i :Integer; x :array [1..N] of Double; y :array [1..N] of Double; d, d1, d2 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 1 to N do begin d := (i - 1) * 1.5 - 4.5; x[i] := d; y[i] := f(d); end; // 0.1刻みで 与えられていない値を補間 for i := 0 to 90 do begin d := i * 0.1 - 4.5; d1 := f(d); d2 := lagrange(d, x, y); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d, d1, d2, d1 - d2])); end; end.
Z:\>fpc Pas0701.pp Free Pascal Compiler version 2.6.2 [2013/02/12] for i386 Copyright (c) 1993-2012 by Florian Klaempfl and others Target OS: Win32 for i386 Compiling Pas0701.pp Linking Pas0701.exe 60 lines compiled, 0.1 sec , 60672 bytes code, 13068 bytes data Z:\>Pas0701 -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
参考文献