ONLY DO WHAT ONLY YOU CAN DO

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

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