Pascalで関数の近似(ネヴィル補間)
与えられた5個の関数値 f(x0), f(x1), f(x2), f(x3), f(x4) を通る4次式は、次のようにして求めることができる.
(1)
(2)
(3)
(4)
(5)
これらの式を使って, 与えられた点以外の点の値を求める.
Program Pas0702(arg); {$MODE delphi} uses SysUtils, Math; const // データ点の数 N = 7; // 元の関数 function f(x:Double):Double; begin result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2); end; // Neville (ネヴィル) 補間 function neville(d:Double; x:array of Double; y:array of Double):Double; var w :array [0..N-1, 0..N-1] of Double; i, j :Integer; begin for i := Low(y) to High(y) do w[0,i] := y[i]; for j := 1 to High(x) do begin for i := 0 to (High(x) - j) do 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; result := w[N-1,0]; end; var i :Integer; x :array [0..N-1] of Double; y :array [0..N-1] of Double; d, d1, d2 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0 to N-1 do begin d := i * 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 := neville(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 Pas0702.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 Pas0702.pp Linking Pas0702.exe 164 lines compiled, 0.1 sec , 60960 bytes code, 13068 bytes data Z:\>Pas0702 -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 以下略