ONLY DO WHAT ONLY YOU CAN DO

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

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