// // interpolation.c // // f(x) = 1 / (1 + x^2) の区間[-1,1]上での等間隔標本点Lagrange補間(反復補間法) // // このファイルをテキストファイルにコピペして,拡張子を.cとして保存してください. // 反復補間法の途中の差分商も表示されます. // f(j, i) = f(x(1), ..., f(x(j-1)), f(x(i)); x) です. // #include #include #define NMAX 100 // double interpolation(double x, double xj[NMAX+1], double fj[NMAX+1], int n_intpol) { double yj[NMAX+1], gj[NMAX+1], ymin, y_tmp, g_tmp; double f_diff[NMAX+1][NMAX+1], fn, eps; int i, j, i1, l_conv; // eps = 1.0e-14; // for (j=1; j<=n_intpol; ++j) { yj[j] = x - xj[j]; gj[j] = fj[j]; } // l_conv = 0; for (i=1; i<=n_intpol; ++i) { ymin = fabs(yj[i]); for (j=i+1; j<=n_intpol; ++j) if (ymin > fabs(yj[j])) { ymin = fabs(yj[j]); i1 = j; } if (i1!=i) { y_tmp = yj[i]; yj[i] = yj[i1]; yj[i1] = y_tmp; g_tmp = gj[i]; gj[i] = gj[i1]; gj[i1] = g_tmp; } // printf("|x(%3d) - x| = %22.15e\n", i, fabs(yj[i])); // f_diff[1][i] = gj[i]; printf("f(%3d, %3d) = %22.15e\n", 1, i, f_diff[1][i]); for (j=2; j<=i; ++j) { f_diff[j][i] = (yj[i] * f_diff[j-1][j-1] - yj[j-1] * f_diff[j-1][i]) / (yj[i] - yj[j-1]); printf("f(%3d, %3d) = %22.15e\n", j, i, f_diff[j][i]); } printf("\n"); if (i>1) if (fabs(f_diff[i][i] - f_diff[i-1][i-1]) < eps) { l_conv = 1; fn = f_diff[i][i]; break; } } // if (l_conv == 0) { printf("Lagrange interpolation not convergent.\n"); fn = f_diff[n_intpol][n_intpol]; } // return fn; } // void data_intpol(double (*func)(double), double xmin, double xmax, int n_intpol, double xj[NMAX+1], double fj[NMAX+1]) { double dx; int j; // dx = (xmax - xmin) / (n_intpol - 1); for (j=1; j<=n_intpol; ++j) { xj[j] = xmin + (j - 1) * dx; fj[j] = func(xj[j]); } // /* for (j=1; j<=n_intpol; ++j) */ /* printf("x(%3d) = %22.15e f(%3d) = %22.15e\n", j, xj[j], j, fj[j]); */ } // double func(double x) { return 1.0 / ( 1.0 + x*x ); } // main() { int n_intpol, k; double xj[NMAX+1], fj[NMAX+1], xmin, xmax; double x, fn, fn_exact, error; double x_array[2] = {.12345,.98765}; // xmin = - 1.0; xmax = 1.0; n_intpol = 55; data_intpol(func, xmin, xmax, n_intpol, xj, fj); // for (k=0; k<=1; ++k) { x = x_array[k]; printf("x = %22.15e\n", x); fn = interpolation(x, xj, fj, n_intpol); fn_exact = func(x); error = fabs(fn - fn_exact); printf("f(x)(interpolation) = %22.15e\n", fn); printf("f(x)(exact) = %22.15e\n", fn_exact); printf("error = %10.3e\n", error); printf("\n"); } // return 0; }