In some cases it makes sense to compute the quadrature formula associated
with the knots. So given x[0], ..., x[n - 1], compute the associated
weights w[0], ..., w[n - 1] and your integral is approximated by w[0] *
f(x[0]) + ... + w[n - 1] * f(x[n - 1]). The snippet below uses GSL to
generate weights so that your quadrature is exact for polynomials of degree
less than n. Check it out, especially if you are using n > 10 since the
system is not well-conditioned. Something to keep in mind is that you
really want nonnegative weights, and there is no guarantee you'll get them.
So if you go this route, make sure you check that none of the generated
weights are negative; if they are, maybe partition the knots into left and
right sets and try again.
Good luck -
- John
#include <stdlib.h>
#include <math.h>
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_vector.h"
static void _getweights( unsigned n, double *x, double *w )
{
unsigned i, j;
double scale;
double t;
double xmax = x[0];
double xmin = x[0];
gsl_matrix *moments = gsl_matrix_alloc( n, n );
gsl_matrix *v = gsl_matrix_alloc( n, n );
gsl_vector *s = gsl_vector_alloc( n );
gsl_vector *work = gsl_vector_alloc( n );
gsl_vector *w0 = gsl_vector_alloc( n );
gsl_vector *r0 = gsl_vector_alloc( n );
for ( i = 1; i < n; i++ ) { /* scale knots to [0, 1] */
if ( x[i] > xmax )
xmax = x[i];
if ( x[i] < xmin )
xmin = x[i];
}
scale = 1 / ( xmax - xmin );
/* Compute the moment matrix and right hand side */
for ( i = 0; i < n; i++ ) {
for ( j = 0; j < n; j++ ) {
t = scale * ( x[j] - xmin );
gsl_matrix_set( moments, i, j, pow( t, ( double ) i ) );
}
gsl_vector_set( r0, i, 1.0 / ( i + 1 ) );
}
gsl_linalg_SV_decomp( moments, v, s, work );
gsl_linalg_SV_solve( moments, v, s, r0, w0 );
for ( i = 0; i < n; i++ )
w[i] = gsl_vector_get( w0, i ) / scale;
gsl_matrix_free( moments );
gsl_matrix_free( v );
gsl_vector_free( s );
gsl_vector_free( work );
gsl_vector_free( w0 );
gsl_vector_free( r0 );
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
