Gauß-Christoffel quadrature

The function gauss_christoffel(N, w, x, total_weight, L, nodes, weights) performs a Gauß-Christoffel quadrature rule to calculate for an arbitrary weight function w(x) (w(x) >= 0, ∀ x) L nodes ( nodes[0] (=x₁), ..., nodes[L-1] (=x_L) ) and weights ( weights[0] (=w₁), ..., weights[L-1] (=w_L) ). The weight function is given to gauss_christoffel in discretized form by the N-dimensional array w with corresponding sampling arguments stored in the array x. The calculated weights conserve the first L moments of the given weight function w(x). The correct scaling of all the moments is ensured by giving to the function the total_weight (=∫ w(x) dx) as an input parameter. For more informations about the Gauß-Christoffel quadrature rule, see https://dlmf.nist.gov/3.5#v.

The Gauß-Christoffel quadrature rule can be used to:

(1) obtain a downsampling of any positive function that needs to conserve the first L moments of the function:

• 0th moment: Σ_i=1,...,L w_i = ∫ w(x) dx,
  
• 1st moment: Σ_i=1,...,L w_i x_i = ∫ w(x)x dx,
  
• ⋮
  
• Lth moment: Σ_i=1,...,L w_i x_i^L = ∫ w(x)x^L dx.

(2) calculate integrals of the form: ∫ w(x) f(x) dx, with w(x) the weight function and f(x) an arbitray function. The nodes x_i and weights w_i obtained from the Gauß-Christoffel quadrature are then used to calculate Σ_i=1,...,L w_i f(x_i) that is a numerical approximation to the integral.

---

Usage

The following example performs a Gauß-Christoffel quadrature for a Gaussian with variance 2 centered around x = 0. An extended version of this example is gauss.c.
Compiling requires linking against a blas and lapack.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "gauss_christoffel.h"


int main(){

    int NN = 1000;                                                          // length of the weight function w
    double sigma2 = 2;                                                      // the variance of the gaussian
    double x0 = 0;                                                          // the mean value of the gaussian
    double x[NN+1], w[NN+1];                                                // allocate arrays for the arguments and the weight function
    for(int k=0; k<=NN; k++){                                               
        x[k] = 5.*((double)k/NN - 0.5);                                     // assign argument values
        w[k] = exp(-(x[k]-x0)*(x[k]-x0)/2./sigma2)/sqrt(2*M_PI*sigma2);     // assign function values for a gaussian
    }

    FILE *F;

    F = fopen("gauss_sample_1000.txt","w+");                                // output the gaussian in the file "gauss_sample_1000.txt"
    for(int k=0; k<=NN; k++)                                                
        fprintf(F, "%f\t%f\n", x[k], w[k]);                                 // write in the first column the arguments and in the second column the values of the weight function
    fclose(F);


    // perform Gauß-Christoffel quadrature with 5 sampling points
    int dim = 5;
    double nodes[dim], weights[dim];

    // call the the Gauß-Christoffel quadrature rule function
    gauss_christoffel(NN+1, w, x, 1., dim, nodes, weights);

    F = fopen("gauss_downsample_N5.txt","w+");                              // output the gaussian in the file "gauss_downsample_N5.txt"
    for(int k=0; k<dim; k++)
        fprintf(F, "%f\t%f\n", nodes[k], weights[k]);                       // write in the first column the nodes and in the second column the weights of the downsampled weight function
    fclose(F);

    return 0;
}