lagrange補間

離散した複数の点を結ぶlagrange補間法。

あんまり実用的でないらしい。

ソース内で複数の点(x,y)を設定し、プログラムを実行することで補間すべき点(x,y)を列挙したファイルlag_dat.i10と補間曲線上の適当な刻み幅の点(x,y)を列挙したファイルlag_dat.o10を作成。

また、それをgnuplot上に描画するためのプロットファイル*.pltを作成する。

lagrange補間は任意の数の点(x,y)を滑らかな曲線でつなぐ関数を求める。

2点なら直線で繋げる、3点なら2次曲線で繋げる、4点なら3次曲線で繋げる、といった事象を利用している。

説明もそこそこに以下リスト。
indent FILENAME -kr -i8
を使うと幸せ。

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

#define N_MAX 3
#define K_MAX 10

double x_fix[N_MAX + 1] = { 0.0, 1.0, 2.0, 3.0 };
double y_fix[N_MAX + 1] = { 1.0, 3.0, 4.0, 2.0 };

char *fname[3] = { "lag_dat.i10", "lag_dat.o10", "lag_dat.plt" };

double p_func(double x_tmp);

int main(void)
{
int k, n;
FILE *fp;
double temp, x, y;

fp = fopen(fname[0], "w");

for (n = 0; n < N_MAX + 1; n++) {
fprintf(fp, "%-lf  ", x_fix[n]);
fprintf(fp, "%-lf\n", y_fix[n]);
}

fclose(fp);

fp = fopen(fname[1], "w");

for (n = 0; n < N_MAX; n++) {
for (k = 0; k < K_MAX; k++) {
temp = (x_fix[n + 1] - x_fix[n]) / (double) K_MAX;
x = x_fix[n] + temp * (double) k;

y = p_func(x);

fprintf(fp, "%-lf  ", x);
fprintf(fp, "%-lf\n", y);
}
}
fprintf(fp, "%-lf  %-lf\n", x_fix[N_MAX], y_fix[N_MAX]);
fclose(fp);

fp = fopen(fname[2], "w");

fprintf(fp,
"plot \"./lag_dat.o10\" using 1:2 with lines,\\\n"
"\"./lag_dat.i10\" using 1:2 with points pt 7\n"
"pause -1\n");
fclose(fp);
}

double p_func(double x_tmp)
{
int m, n;
double temp;
double q[N_MAX][N_MAX + 1];

for (n = 0; n < N_MAX; n++) {
if (n == 0) {
for (m = n + 1; m < N_MAX + 1; m++) {
temp = y_fix[n] * (x_fix[m] - x_tmp)
   - y_fix[m] * (x_fix[n] - x_tmp);
q[n][m] = temp / (x_fix[m] - x_fix[n]);
}
} else {
for (m = n + 1; m < N_MAX + 1; m++) {
temp = q[n - 1][n] * (x_fix[m] - x_tmp)
   - q[n - 1][m] * (x_fix[n] - x_tmp);
q[n][m] = temp / (x_fix[m] - x_fix[n]);
}
}
}
return q[N_MAX - 1][N_MAX];
}