ポリゴン関数による多項式補間法のC言語実装

参考文献

数値計算の基礎と応用[新訂版]
数値解析学への入門
杉浦　洋(南山大学教授)　著
発行日 2009/12/10

準備

オンラインコンパイラを使用します。

ソースコード

sample.c

#include <stdio.h>
#include <math.h>

#define M 7

void NewtonCoef(double xi[], int m, double b[]) {
    for (int n = 0; n <= m; n++) {
        b[n] = exp(xi[n]);
        for (int l = 0; l < n; l++) {
            b[n] = (b[n] - b[l]) / (xi[n] - xi[l]);
        }
    }
}

double p(double x, int n, double xi[], double b[]) {
    double y = b[n];
    for (int l = n - 1; l >= 0; l--) {
        y = (x - xi[l]) * y + b[l];
    }
    return y;
}

double f(double x) {
    return exp(x);
}

int main() {
    int m = M;
    double Pi = 3.14159265358979323846;
    double dt = Pi / (m + 1);
    double xi[M + 1];
    double b[M + 1];

    for (int i = 0; i <= m; i++) {
        xi[i] = 0.5 * pow((i + 0.5) * dt, 2);
    }

    NewtonCoef(xi, m, b);

    printf("degree=%d\n", m);
    int npoints = 10;
    double dx = 1.0 / npoints;

    for (int i = 0; i <= npoints; i++) {
        double x = -0.5 + i * dx;
        double y = p(x, m, xi, b);
        printf("p(%4.1f)=%17.10e error=%9.2e\n", x, y, y - f(x));
    }

    return 0;
}

実行結果

console

degree=7
p(-0.5)= 6.0053261664e-01 error=-6.00e-03
p(-0.4)= 6.6727571513e-01 error=-3.04e-03
p(-0.3)= 7.3943815181e-01 error=-1.38e-03
p(-0.2)= 8.1820681408e-01 error=-5.24e-04
p(-0.1)= 9.0469626135e-01 error=-1.41e-04
p( 0.0)= 9.9999065019e-01 error=-9.35e-06
p( 0.1)= 1.1051809831e+00 error= 1.01e-05
p( 0.2)= 1.2213982646e+00 error=-4.49e-06
p( 0.3)= 1.3498431463e+00 error=-1.57e-05
p( 0.4)= 1.4918126465e+00 error=-1.21e-05
p( 0.5)= 1.6487245249e+00 error= 3.25e-06