補間ライブラリ 使用例



[ 簡単な説明 ]

補間ライブラリの使用例です。

最小2乗近似 lstsq( ) では、あらかじめ設定した3次の多項式から得た8組の標本点から、3次多項式の係数を計算しています。

ラグランジュ補間 lagra( ) では、指数関数 exp( ) から得られた5組の標本点から、補間値を計算し、真値と比較しています。

スプライン補間 splint( ) では、sin 関数から得られた20組の標本点から、補間値を計算し、真値と比較しています。

チェビシェフ近似 chebyshev( ) では、指数関数 exp( )を範囲 -2 〜 2 で3次、7次、15次の多項式で近似し、真値と比較しています。次数が高くなると、級数展開に近い係数となっているのが判ります。

プログラム・ソース("test18.c")           top (先頭に戻る)
/*		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;
}

出力結果           top (先頭に戻る)
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