[ 簡単な説明 ] 補間ライブラリの使用例です。 最小2乗近似 lstsq( ) では、あらかじめ設定した3次の多項式から得た8組の標本点から、3次多項式の係数を計算しています。 ラグランジュ補間 lagra( ) では、指数関数 exp( ) から得られた5組の標本点から、補間値を計算し、真値と比較しています。 スプライン補間 splint( ) では、sin 関数から得られた20組の標本点から、補間値を計算し、真値と比較しています。 チェビシェフ近似 chebyshev( ) では、指数関数 exp( )を範囲 -2 〜 2 で3次、7次、15次の多項式で近似し、真値と比較しています。次数が高くなると、級数展開に近い係数となっているのが判ります。 |
/* test18.c */ #include <stdio.h> #include "sslib.h" double _f(double x); void print_coef(double a[], int n); void eval(double a[], int n, double min, double max, int m); double _f(double x) { return exp(x); } void print_coef(double a[], int n) { int i; printf("f(x) = %9.6f * X^%d\n", a[n], n); for(i = n - 1; i >= 2; i--) { printf((a[i] >= 0.)? " + ": " "); printf("%9.6f * X^%d\n", a[i], i); } printf((a[1] >= 0.)? " + ": " "); printf("%9.6f * X\n", a[1]); printf((a[0] >= 0.)? " + ": " "); printf("%9.6f\n\n", a[0]); } void eval(double a[], int n, double min, double max, int m) { double h, w, x, y; int i; printf("*** Evaluation ***\n"); h = 1. / (double)m; for(i = 0; i <= m; i++) { x = min + (max - min) * h * i; y = _f(x); w = polynomial(a, x, n); printf("x=%4.1f f(x)=%8.6f approximation=%8.6f error=%11.3e\n", x, y, w, w - y); } printf("\n"); } int main(void) { static double x[50] = { -3., -2., -1., 0., 1., 2., 3., 4.}; static double y[50] = { -58., -26., -10., -4., -2., 2., 14., 40.}; static double x1[5] = { 0.15, 0.20, 0.25, 0.30, 0.35 }; static double y1[5] = { 0.86070798, 0.81873075, 0.77880078, 0.74081822, 0.70468809}; double a[20], c[4]; double max, min, w, xi, yi, z; int i, m, n; n = 8; /* データ個数 */ m = 3; /* 近似多項式次数 */ printf("y = x ^ 3 - 2 * x ^ 2 + 3 * x - 4\n"); lstsq(x, y, n, m, c); printf("y = (%6.3f) * x ^ 3 + (%6.3f) * x ^ 2 + (%6.3f) * x + (%6.3f)\n", c[3], c[2], c[1], c[0]); n = 5; xi = 0.27; printf(" x lagra exp\n"); z = exp(-xi); w = lagra(x1, y1, n, xi); printf("%4.2f %14.7e %14.7e\n", xi, w, z); n = 20; for(i = 0; i <= n; i++) { x[i] = M_PI * i / n; y[i] = sin(x[i]); } xi = 1.0; yi = splint(x, y, n, xi); z = sin(xi); printf(" Interpolation by Cubic Spline function\n"); printf(" Interpolation = %10.6f at x = %8.3f\n", yi, xi); printf(" Exact = %10.6f\n", z); printf(" Error = %10.6f\n", yi - z); printf("\n*** Chebyshev Approximation ***\n\n"); max = 2; min = -2; n = 3; chebyshev(min, max, n, a); print_coef(a, n); eval(a, n, min, max, 10); n = 7; chebyshev(min, max, n, a); print_coef(a, n); eval(a, n, min, max, 10); n = 15; chebyshev(min, max, n, a); print_coef(a, n); eval(a, n, min, max, 10); return 1; } |
y = x ^ 3 - 2 * x ^ 2 + 3 * x - 4 y = ( 1.000) * x ^ 3 + (-2.000) * x ^ 2 + ( 3.000) * x + (-4.000) x lagra exp 0.27 7.6337949e-01 7.6337949e-01 Interpolation by Cubic Spline function Interpolation = -1.102539427130981640000000000000000000000e+214 at x = 1.000 Exact = 0.841471 Error = -1.102539427130981640000000000000000000000e+214 *** Chebyshev Approximation *** f(x) = 0.202914 * X^3 + 0.687348 * X^2 + 0.981666 * X + 0.904834 *** Evaluation *** x=-2.0 f(x)=0.135335 approximation=0.067579 error= -6.776e-02 x=-1.6 f(x)=0.201897 approximation=0.262642 error= 6.075e-02 x=-1.2 f(x)=0.301194 approximation=0.365980 error= 6.479e-02 x=-0.8 f(x)=0.449329 approximation=0.455511 error= 6.182e-03 x=-0.4 f(x)=0.670320 approximation=0.609157 error= -6.116e-02 x= 0.0 f(x)=1.000000 approximation=0.904834 error= -9.517e-02 x= 0.4 f(x)=1.491825 approximation=1.420463 error= -7.136e-02 x= 0.8 f(x)=2.225541 approximation=2.233962 error= 8.421e-03 x= 1.2 f(x)=3.320117 approximation=3.423251 error= 1.031e-01 x= 1.6 f(x)=4.953032 approximation=5.066248 error= 1.132e-01 x= 2.0 f(x)=7.389056 approximation=7.240873 error= -1.482e-01 f(x) = 0.000222 * X^7 + 0.001600 * X^6 + 0.008274 * X^5 + 0.041129 * X^4 + 0.166714 * X^3 + 0.500433 * X^2 + 0.999994 * X + 0.999946 *** Evaluation *** x=-2.0 f(x)=0.135335 approximation=0.135291 error= -4.433e-05 x=-1.6 f(x)=0.201897 approximation=0.201877 error= -1.942e-05 x=-1.2 f(x)=0.301194 approximation=0.301174 error= -2.019e-05 x=-0.8 f(x)=0.449329 approximation=0.449378 error= 4.925e-05 x=-0.4 f(x)=0.670320 approximation=0.670322 error= 2.080e-06 x= 0.0 f(x)=1.000000 approximation=0.999946 error= -5.420e-05 x= 0.4 f(x)=1.491825 approximation=1.491827 error= 2.270e-06 x= 0.8 f(x)=2.225541 approximation=2.225600 error= 5.867e-05 x= 1.2 f(x)=3.320117 approximation=3.320091 error= -2.627e-05 x= 1.6 f(x)=4.953032 approximation=4.953005 error= -2.761e-05 x= 2.0 f(x)=7.389056 approximation=7.388987 error= -6.890e-05 f(x) = 0.000000 * X^15 + 0.000000 * X^14 + 0.000000 * X^13 + 0.000000 * X^12 + 0.000000 * X^11 + 0.000000 * X^10 + 0.000003 * X^9 + 0.000025 * X^8 + 0.000198 * X^7 + 0.001389 * X^6 + 0.008333 * X^5 + 0.041667 * X^4 + 0.166667 * X^3 + 0.500000 * X^2 + 1.000000 * X + 1.000000 *** Evaluation *** x=-2.0 f(x)=0.135335 approximation=0.135335 error= 9.361e-13 x=-1.6 f(x)=0.201897 approximation=0.201897 error= -3.593e-13 x=-1.2 f(x)=0.301194 approximation=0.301194 error= 2.463e-12 x=-0.8 f(x)=0.449329 approximation=0.449329 error= -3.085e-12 x=-0.4 f(x)=0.670320 approximation=0.670320 error= 3.401e-12 x= 0.0 f(x)=1.000000 approximation=1.000000 error= -3.559e-12 x= 0.4 f(x)=1.491825 approximation=1.491825 error= 3.507e-12 x= 0.8 f(x)=2.225541 approximation=2.225541 error= -3.559e-12 x= 1.2 f(x)=3.320117 approximation=3.320117 error= 4.687e-12 x= 1.6 f(x)=4.953032 approximation=4.953032 error= -6.449e-12 x= 2.0 f(x)=7.389056 approximation=7.389056 error= 6.495e-12 |