ベッセル関数による多項式補間法のC言語実装

参考ページ

準備

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

ソースコード

sample.c

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

#define M 7
#define NPOINTS 10
#define PI 3.14159265358979323846

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

double bessel_j0(double x, int num_terms) {
    double result = 0.0;
    for (int k = 0; k < num_terms; k++) {
        double term = pow(-1, k) / (tgamma(k + 1) * tgamma(k + 1)) * pow(x / 2, 2 * k);
        result += term;
    }
    return result;
}

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

    // ベッセル関数を使用
    for (int i = 0; i <= m; i++) {
        xi[i] = bessel_j0(i * dt, 50);
    }

    NewtonCoef(xi, m, b);
    printf("degree=%d\n", m);

    double dx = 1.0 / NPOINTS;
    double x_vals[NPOINTS + 1];
    double y_vals[NPOINTS + 1];
    double y_actual[NPOINTS + 1];

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

    return 0;
}

実行結果

console

degree=7
p(-0.5)= 6.0651391276e-01 error=-1.67e-05
p(-0.4)= 6.7031453552e-01 error=-5.51e-06
p(-0.3)= 7.4081690019e-01 error=-1.32e-06
p(-0.2)= 8.1873063103e-01 error=-1.22e-07
p(-0.1)= 9.0483747139e-01 error= 5.34e-08
p( 0.0)= 1.0000000106e+00 error= 1.06e-08
p( 0.1)= 1.1051709038e+00 error=-1.43e-08
p( 0.2)= 1.2214027522e+00 error=-5.97e-09
p( 0.3)= 1.3498588124e+00 error= 4.78e-09
p( 0.4)= 1.4918247018e+00 error= 4.16e-09
p( 0.5)= 1.6487212693e+00 error=-1.38e-09

[Execution complete with exit code 0]