C++ - ニュートン補間!

Updated:


何組かの x, y データが与えられ、これらの点全てを通る補間多項式を求める方法に「ニュートン補間」というものがあります。

先日は「ラグランジュ補間」について紹介しました。

以下、一部 \(\TeX\) で記載。

INTERPOLATE_NEWTON_1

INTERPOLATE_NEWTON_2

アルゴリズムとしては、係数(差分商)を求め、求まった多項式から補間点を算出していく形になる。

以下、C++ によるサンプルソースです。

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2
  • ニュートン補間そのものについての詳細は割愛。

1. C++ ソース作成

File: interpolate_newton.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
/*********************************************
 * ニュートン補間
 *********************************************/
#include <iostream>  // for cout
#include <stdio.h>   // for printf()

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 各種変数
    int    n;       // あらかじめ与える点の数
    int    i, j;    // LOOP インデックス
    double t;       // LOOP インデックス
    double c[100];  // 係数配列
    double w[100];  // 作業用配列
    double s;       // 総和

    public:
        // 計算
        void calc();
        // ニュートン補間
        double interpolateNewton(double x[], double y[], int n, double t);
};

/*
 * 計算
 */
void Calc::calc()
{
    // あらかじめ与える点
    double x[] = {0.0, 2.0, 3.0, 5.0, 8.0};
    double y[] = {0.8, 3.2, 2.8, 4.5, 1.9};

    // 点の数
    n = sizeof(x) / sizeof(x[0]);

    // 結果出力
    printf("      x      y\n");
    for (t = x[0]; t <= x[n - 1]; t += .5)
        printf("%7.2f%7.2f\n", t, interpolateNewton(x, y, n, t));
}

/*
 * ニュートン補間
 */
double Calc::interpolateNewton(double x[], double y[], int n, double t)
{
    // 係数
    for (i = 0; i < n; i++) {
        w[i] = y[i];
        for (j = i - 1; j >= 0; j--)
            w[j] = (w[j+1] - w[j]) / (x[i] - x[j]);
        c[i] = w[0];
    }

    // 総和
    s = c[n-1];
    for (i = n - 2; i >= 0; i--)
        s = s * (t - x[i]) + c[i];

    return s;
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // ニュートン補間計算
        objCalc.calc();
    }
    catch (...) {
        cout << "例外発生!" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

2. C++ ソースコンパイル

$ g++ interpolate_newton.cpp -o interpolate_newton

何も出力されなければ成功です。

3. 実行

$ ./interpolate_newton
      x      y
   0.00   0.80
   0.50   2.49
   1.00   3.23
   1.50   3.37
   2.00   3.20
   2.50   2.95
   3.00   2.80
   3.50   2.85
   4.00   3.17
   4.50   3.74
   5.00   4.50
   5.50   5.32
   6.00   6.03
   6.50   6.37
   7.00   6.05
   7.50   4.70
   8.00   1.90

あらかじめ与えられた点以外の点も計算されているのが分かる。
また、「ラグランジュ補間」と同じ結果になっているのも分かる。(C++ - ラグランジュ補間!

4. グラフ

参考までに、 R でグラフを作成してみた。

INTERPOLATE_GRAPH


数学はおもしろいけど、コンピュータで実証するというのもおもしろいです。

※ちなみに最近の当方の C++ アルゴリズムについての記事は、古い C によるアルゴリズムに関する書物を参考に C++ に移植した形態となっています。

以上。





 

Sponsored Link

 

Comments