タンジェントによる多項式補間法のC言語実装

参考文献

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

準備

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

ソースコード

sample.c

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

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);
}

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

int main() {
    int m = 7;
    double Pi = M_PI;
    double dt = Pi / (m + 1);
    double xi[m + 1];
    double b[m + 1];
    
    for (int i = 0; i <= m; i++) {
        xi[i] = 0.5 * tan((i + 0.5) * dt);
    }
    
    NewtonCoef(xi, m, b);
    printf("degree=%d\n", m);
    
    int npoints = 100;
    double dx = 1.0 / npoints;
    double x, y, y_actual;
    
    for (int i = 0; i <= npoints; i++) {
        x = -0.5 + i * dx;
        y = p(x, m, xi, b);
        y_actual = f(x);
        printf("p(%4.1f)=%17.10e error=%9.2e\n", x, y, y - exp(x));
    }
    
    return 0;
}

実行結果

console

degree=7
p(-0.5)= 6.0652907050e-01 error=-1.59e-06
p(-0.5)= 6.1262493040e-01 error=-1.46e-06
p(-0.5)= 6.1878205271e-01 error=-1.34e-06
p(-0.5)= 6.2500105220e-01 error=-1.22e-06
p(-0.5)= 6.3128254996e-01 error=-1.10e-06
p(-0.5)= 6.3762717339e-01 error=-9.78e-07
p(-0.4)= 6.4403555630e-01 error=-8.65e-07
p(-0.4)= 6.5050833896e-01 error=-7.56e-07
p(-0.4)= 6.5704616817e-01 error=-6.52e-07
p(-0.4)= 6.6364969729e-01 error=-5.53e-07
p(-0.4)= 6.7031958633e-01 error=-4.60e-07
p(-0.4)= 6.7705650203e-01 error=-3.72e-07
p(-0.4)= 6.8386111786e-01 error=-2.91e-07
p(-0.4)= 6.9073411415e-01 error=-2.16e-07
p(-0.4)= 6.9767617810e-01 error=-1.48e-07
p(-0.3)= 7.0468800390e-01 error=-8.58e-08
p(-0.3)= 7.1177029274e-01 error=-3.00e-08
p(-0.3)= 7.1892375293e-01 error= 1.95e-08
p(-0.3)= 7.2614909991e-01 error= 6.28e-08
p(-0.3)= 7.3344705638e-01 error= 1.00e-07
p(-0.3)= 7.4081835232e-01 error= 1.32e-07
p(-0.3)= 7.4826372507e-01 error= 1.57e-07
p(-0.3)= 7.5578391943e-01 error= 1.78e-07
p(-0.3)= 7.6337968768e-01 error= 1.93e-07
p(-0.3)= 7.7105178970e-01 error= 2.04e-07
p(-0.2)= 7.7880099303e-01 error= 2.10e-07
p(-0.2)= 7.8662807292e-01 error= 2.12e-07
p(-0.2)= 7.9453381243e-01 error= 2.10e-07
p(-0.2)= 8.0251900249e-01 error= 2.05e-07
p(-0.2)= 8.1058444200e-01 error= 1.96e-07
p(-0.2)= 8.1873093786e-01 error= 1.85e-07
p(-0.2)= 8.2695930511e-01 error= 1.71e-07
p(-0.2)= 8.3527036697e-01 error= 1.56e-07
p(-0.2)= 8.4366495490e-01 error= 1.38e-07
p(-0.2)= 8.5214390875e-01 error= 1.20e-07
p(-0.1)= 8.6070807677e-01 error= 1.00e-07
p(-0.1)= 8.6935831574e-01 error= 8.03e-08
p(-0.1)= 8.7809549101e-01 error= 6.01e-08
p(-0.1)= 8.8692047663e-01 error= 3.99e-08
p(-0.1)= 8.9583415542e-01 error= 2.01e-08
p(-0.1)= 9.0483741905e-01 error= 1.01e-09
p(-0.1)= 9.1393116811e-01 error=-1.72e-08
p(-0.1)= 9.2311631225e-01 error=-3.41e-08
p(-0.1)= 9.3239377021e-01 error=-4.97e-08
p(-0.1)= 9.4176446995e-01 error=-6.36e-08
p(-0.0)= 9.5122934873e-01 error=-7.58e-08
p(-0.0)= 9.6078935321e-01 error=-8.59e-08
p(-0.0)= 9.7044543953e-01 error=-9.40e-08
p(-0.0)= 9.8019857341e-01 error=-9.99e-08
p(-0.0)= 9.9004973024e-01 error=-1.04e-07
p( 0.0)= 9.9999989521e-01 error=-1.05e-07
p( 0.0)= 1.0100500633e+00 error=-1.04e-07
p( 0.0)= 1.0202012397e+00 error=-1.00e-07
p( 0.0)= 1.0304544393e+00 error=-9.46e-08
p( 0.0)= 1.0408106875e+00 error=-8.67e-08
p( 0.1)= 1.0512710198e+00 error=-7.66e-08
p( 0.1)= 1.0618364821e+00 error=-6.45e-08
p( 0.1)= 1.0725081308e+00 error=-5.05e-08
p( 0.1)= 1.0832870329e+00 error=-3.47e-08
p( 0.1)= 1.0941742662e+00 error=-1.75e-08
p( 0.1)= 1.1051709191e+00 error= 1.04e-09
p( 0.1)= 1.1162780911e+00 error= 2.06e-08
p( 0.1)= 1.1274968926e+00 error= 4.10e-08
p( 0.1)= 1.1388284451e+00 error= 6.18e-08
p( 0.1)= 1.1502738817e+00 error= 8.28e-08
p( 0.2)= 1.1618343464e+00 error= 1.04e-07
p( 0.2)= 1.1735109951e+00 error= 1.24e-07
p( 0.2)= 1.1853049949e+00 error= 1.44e-07
p( 0.2)= 1.1972175249e+00 error= 1.62e-07
p( 0.2)= 1.2092497761e+00 error= 1.78e-07
p( 0.2)= 1.2214029512e+00 error= 1.93e-07
p( 0.2)= 1.2336782652e+00 error= 2.05e-07
p( 0.2)= 1.2460769452e+00 error= 2.15e-07
p( 0.2)= 1.2586002307e+00 error= 2.21e-07
p( 0.2)= 1.2712493736e+00 error= 2.23e-07
p( 0.2)= 1.2840256385e+00 error= 2.22e-07
p( 0.3)= 1.2969303025e+00 error= 2.16e-07
p( 0.3)= 1.3099646559e+00 error= 2.05e-07
p( 0.3)= 1.3231300016e+00 error= 1.89e-07
p( 0.3)= 1.3364276558e+00 error= 1.68e-07
p( 0.3)= 1.3498589482e+00 error= 1.41e-07
p( 0.3)= 1.3634252213e+00 error= 1.07e-07
p( 0.3)= 1.3771278317e+00 error= 6.74e-08
p( 0.3)= 1.3909681494e+00 error= 2.10e-08
p( 0.3)= 1.4049475582e+00 error=-3.23e-08
p( 0.3)= 1.4190674559e+00 error=-9.27e-08
p( 0.4)= 1.4333292544e+00 error=-1.60e-07
p( 0.4)= 1.4477343799e+00 error=-2.35e-07
p( 0.4)= 1.4622842728e+00 error=-3.17e-07
p( 0.4)= 1.4769803882e+00 error=-4.06e-07
p( 0.4)= 1.4918241958e+00 error=-5.02e-07
p( 0.4)= 1.5068171803e+00 error=-6.05e-07
p( 0.4)= 1.5219608411e+00 error=-7.14e-07
p( 0.4)= 1.5372566931e+00 error=-8.30e-07
p( 0.4)= 1.5527062662e+00 error=-9.52e-07
p( 0.5)= 1.5683111059e+00 error=-1.08e-06
p( 0.5)= 1.5840727732e+00 error=-1.21e-06
p( 0.5)= 1.5999928452e+00 error=-1.35e-06
p( 0.5)= 1.6160729145e+00 error=-1.49e-06
p( 0.5)= 1.6323145902e+00 error=-1.63e-06
p( 0.5)= 1.6487194974e+00 error=-1.77e-06

[Execution complete with exit code 0]