C++で関数の近似(スプライン補間)

をスプライン補間で近似する

複数の点が与えられているとき, これらを１つの多項式でつなげるのではなく,
２点ずつの区間に分けて, 小区間ごとに別々の３次式を設定する.

(1) 連立1次方程式をたてる

(2) 連立1次方程式を解いて s_i を求める

(3) 区間 x_i 〜 x_i+1 の３次式は次のようになる.

この式を使って, 与えられた点以外の点の値を求める.

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

// データ点の数
const int  N = 7; 

// 元の関数
double f(double x);

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[]);

int main()
{
    double x[N], y[N];

    // 1.5刻みで -4.5〜4.5 まで, ７点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i] = d;
        y[i] = f(d);
    }

    // ３項方程式の係数の表を作る
    double a[N], b[N], c[N], d[N];
    for (int i = 1; i < N - 1; i++)
    {
        a[i] =         x[i]   - x[i-1]; 
        b[i] = 2.0 *  (x[i+1] - x[i-1]); 
        c[i] =         x[i+1] - x[i]; 
        d[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1]));
    }
    // ３項方程式を解く (ト−マス法)
    double g[N], s[N];
    g[1] = b[1];
    s[1] = d[1];
    for (int i = 2; i < N - 1; i++)
    {
        g[i] = b[i] - a[i] * c[i-1] / g[i-1];
        s[i] = d[i] - a[i] * s[i-1] / g[i-1];
    }
    double z[N];
    z[N-2] = s[N-2] / g[N-2];
    for (int i = N - 3; i >= 1; i--)
    {
        z[i] = (s[i] - c[i] * z[i+1]) / g[i];
    }

    // 0.1刻みで 与えられていない値を補間 
    for (int i = 0; i <= 80; i++)
    {
        double d  = i * 0.1 - 4.0;
        double d1 = spline(d, x, y, z);
        double d2 = f(d);

        // 元の関数と比較
        cout << setw(5) << fixed << setprecision(2) << d       << " : "; 
        cout << setw(8) << fixed << setprecision(5) << d1      << " - "; 
        cout << setw(8) << fixed << setprecision(5) << d2      << " = "; 
        cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl; 
    }
   return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3)/6 + pow(x,5)/120;
}

// Spline (スプライン) 補間
double spline(double d, double x[], double y[], double z[])
{
    // 補間関数値がどの区間にあるか
    int k=-1;
    for (int i = 0; i < N - 1; i++)
    {
        if (x[i] <= d && d < x[i+1])
        {
            k = i;
            break;
        }
    }
    if (k < 0) return -1;

    double d1 = x[k+1] - d;
    double d2 = d      - x[k];
    double d3 = x[k+1] - x[k];
    return      (z[k] * pow(d1,3) + z[k+1] * pow(d2,3)) / (6.0 * d3)
              + (y[k]   / d3 - z[k]   * d3 / 6.0) * d1  
              + (y[k+1] / d3 - z[k+1] * d3 / 6.0) * d2;

}

Z:\>bcc32 CP0705.cpp
Borland C++ 5.5.1 for Win32 Copyright (c) 1993, 2000 Borland
cp0705.cpp:
Turbo Incremental Link 5.00 Copyright (c) 1997, 2000 Borland

Z:\>CP0705
-4.00 : -2.90573 - -1.86667 = -1.03906
-3.90 : -2.57503 - -1.53218 = -1.04285
-3.80 : -2.25858 - -1.25760 = -1.00099
-3.70 : -1.95876 - -1.03650 = -0.92226
-3.60 : -1.67794 - -0.86285 = -0.81509
-3.50 : -1.41849 - -0.73099 = -0.68750
-3.40 : -1.18279 - -0.63562 = -0.54717
-3.30 : -0.97322 - -0.57178 = -0.40144
-3.20 : -0.79215 - -0.53487 = -0.25728
-3.10 : -0.64195 - -0.52060 = -0.12135
-3.00 : -0.52500 - -0.52500 =  0.00000
-2.90 : -0.44268 - -0.54443 =  0.10175
-2.80 : -0.39235 - -0.57553 =  0.18318
-2.70 : -0.37041 - -0.61524 =  0.24484
-2.60 : -0.37321 - -0.66078 =  0.28757
-2.50 : -0.39714 - -0.70964 =  0.31250
-2.40 : -0.43856 - -0.75955 =  0.32099
-2.30 : -0.49386 - -0.80853 =  0.31466
-2.20 : -0.55942 - -0.85480 =  0.29539
以下略

ONLY DO WHAT ONLY YOU CAN DO

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

参考文献