微積分ライブラリ 使用例

/*		test4.c		*/
#include	<stdio.h>
#include	"sslib.h"

/* should be defind by user as _f(x) for cheb3(), cheb4(), chev6()
									 for dgl3(), dgl10(), dgl20(), dgl32(), dgl48()	*/

double _f(double x)
{
	return exp(x);
}

int main(void)
{
	int i, n;
	double x[100], y[100], a, b, exa, h, *p, *q, s;

	n = 100;
	for(i = 0, p = x, q = y; i < n; i++, p++, q++)
	{
		*p = (double)i / (double)(n - 1);
		*q = exp(*p);
	}
	s = trap(x, y, n);
	exa = exp(1.) - 1.;
	printf("* Integral by Trapezoid's rule\n");
	printf("    exp(x)    x = 0 --> 1  (devided point %d)\n", n);
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	n = 21;
	h = M_PI / (double)(n - 1);
	for(i = 0, p = y; i < n; i++, p++)	*p = sin(h * (double)i);
	s = simpei(y, n, h);
	exa = 2.;
	printf("* Integral by Simpson's method\n");
	printf("    exp(x)    x = 0 --> 1  (devided point %d)\n", n);
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	n = 21;
	h = M_PI / (double)(n - 1);
	for(i = 0, p = x, q = y; i < n; i++, p++, q++)
	{
		*p = h * (double)i;
		if ((i != 0) && (i != n - 1))	*p += .005 * (double)i;
		*q = sin(*p);
	}
	s = simpui(x, y, n);
	exa = 2.;
	printf("* Integral by Simpson's method (unequal interval)\n");
	printf("    exp(x)    x = 0 --> 1  (devided point %d)\n", n);
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	n = 21;
	h = M_PI / (double)(n - 1);
	for(i = 0, p = x, q = y; i < n; i++, p++, q++)
	{
		*p = h * (double)i;
		*q = sin(*p);
	}
	s = splitg(x, y, n);
	exa = 2.;
	printf("* Integral by Cubic Spline Function's rule\n");
	printf("    exp(x)    x = 0 --> 1  (devided point %d)\n", n);
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	a = 0.;
	b = 1.;
	printf("* Integration by Chebyshev's method\n");
	printf("   exp(x)   x=%4.1f --> %4.1f\n", a, b);
	exa = exp(1.) - 1.;
	s = cheb3(a, b);
	printf("*   3 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = cheb4(a, b);
	printf("*   4 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = cheb6(a, b);
	printf("*   6 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	a = 0.;
	b = 1.;
	printf("* Integration by Gauss-Legendre's method\n");
	printf("   exp(x)   x=%4.1f --> %4.1f\n", a, b);
	exa = exp(1.) - 1.;
	s = dgl3(a, b);
	printf("*   3 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = dgl10(a, b);
	printf("*   10 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = dgl20(a, b);
	printf("*   20 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = dgl32(a, b);
	printf("*   32 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);
	s = dgl48(a, b);
	printf("*   48 - point\n");
	printf("    Integration = %22.15e\n", s);
	printf("    Exact       = %22.15e\n", exa);
	printf("    Error       = %12.5e\n", s - exa);

	return 1;
}