Blame doc/examples/integration2.c
|
Packit |
67cb25 |
#include <stdio.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_integration.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
f(double x, void * params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int m = *(int *) params;
|
|
Packit |
67cb25 |
double f = gsl_pow_int(x, m) + 1.0;
|
|
Packit |
67cb25 |
return f;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main (int argc, char *argv[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_integration_fixed_workspace * w;
|
|
Packit |
67cb25 |
const gsl_integration_fixed_type * T = gsl_integration_fixed_hermite;
|
|
Packit |
67cb25 |
int m = 10;
|
|
Packit |
67cb25 |
int n = 6;
|
|
Packit |
67cb25 |
double expected, result;
|
|
Packit |
67cb25 |
gsl_function F;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (argc > 1)
|
|
Packit |
67cb25 |
m = atoi(argv[1]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (argc > 2)
|
|
Packit |
67cb25 |
n = atoi(argv[2]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = gsl_integration_fixed_alloc(T, n, 0.0, 1.0, 0.0, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F.function = &f;
|
|
Packit |
67cb25 |
F.params = &m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_integration_fixed(&F, &result, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (m % 2 == 0)
|
|
Packit |
67cb25 |
expected = M_SQRTPI + gsl_sf_gamma(0.5*(1.0 + m));
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
expected = M_SQRTPI;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
printf ("m = %d\n", m);
|
|
Packit |
67cb25 |
printf ("intervals = %zu\n", gsl_integration_fixed_n(w));
|
|
Packit |
67cb25 |
printf ("result = % .18f\n", result);
|
|
Packit |
67cb25 |
printf ("exact result = % .18f\n", expected);
|
|
Packit |
67cb25 |
printf ("actual error = % .18f\n", result - expected);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_integration_fixed_free (w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|