ONLY DO WHAT ONLY YOU CAN DO

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

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.54.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
以下略
参考文献