ニュートン補間法によるチェビシェフ補間のC言語実装

参考文献

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

参考ページ

準備

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

ソースコード

sample.c

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

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

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

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

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


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

実行結果

console

degree=7
p(-0.5)= 6.0653065899e-01 error=-7.21e-10
p(-0.4)= 6.7032004573e-01 error=-3.07e-10
p(-0.3)= 7.4081822037e-01 error=-3.11e-10
p(-0.2)= 8.1873075381e-01 error= 7.36e-10
p(-0.1)= 9.0483741807e-01 error= 3.02e-11
p( 0.0)= 9.9999999924e-01 error=-7.61e-10
p( 0.1)= 1.1051709181e+00 error= 3.08e-11
p( 0.2)= 1.2214027589e+00 error= 7.70e-10
p( 0.3)= 1.3498588072e+00 error=-3.32e-10
p( 0.4)= 1.4918246973e+00 error=-3.36e-10
p( 0.5)= 1.6487212699e+00 error=-8.06e-10