Pascalで関数の近似(ニュートン補間)
をニュートン補間で近似する
与えられた4個の関数値 f(x0), f(x1), f(x2), f(x3) を通る3次式を求める場合, まず次のような表を作る.
このとき,
と定義する.
これを x0とx1 の第1差分商, または1階差分商という.
同様に
と定義し, これを x0, x1, x2 の第2差分商, または2階差分商という.
第n差分商を
で表すと, 与えられた n+1点を通る n次式は次のように表すことができる.
この式を使って, 与えられた点以外の点の値を求める.
Program Pas0703(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; // Newton (ニュートン) 補間 function newton(d:Double; x:array of Double; a:array of Double):Double; var sum, prod :Double; i, j :Integer; begin sum := a[0]; for i := 1 to High(x) do begin prod := 1; for j := Low(x) to i-1 do prod := prod * (d - x[j]); sum := sum + (a[i] * prod); end; result := sum; end; var i, j :Integer; x :array [0..N-1] of Double; y :array [0..N-1] of Double; a :array [0..N-1] of Double; d :array [0..N-1, 0..N-1] of Double; d3, d1, d2 :Double; begin // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット for i := 0 to N-1 do begin d3 := i * 1.5 - 4.5; x[i] := d3; y[i] := f(d3); end; // 差分商の表を作る for j := 0 to N-1 do d[0,j] := y[j]; for i := 1 to N-1 do begin for j := 0 to N-i-1 do d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (x[j+i] - x[j]); end; // n階差分商 for j := 0 to N-1 do a[j] := d[j,0]; // 0.1刻みで 与えられていない値を補間 for i := 0 to 90 do begin d3 := i * 0.1 - 4.5; d1 := f(d3); d2 := newton(d3, x, a); // 元の関数と比較 writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d3, d1, d2, d1 - d2])); end; end.
Z:\>fpc Pas0703.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 Pas0703.pp Linking Pas0703.exe 87 lines compiled, 0.2 sec , 61024 bytes code, 13068 bytes data Z:\>Pas0703 -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 以下略