|
[ 簡単な説明 ] 補間ライブラリの使用例です。 最小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
|