Re: Numerical quadrature performed on a set of points

[crow@.......]

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 );
}