Update: Gauss Numeric Integrator as open source project hosted on codeplex: http://gauss.codeplex.com

In this post it will be presented C# implementation of numerical calculation of  integral within interval based on Gauss Quadrature method.

As you know definite integral can be expressed as:

.

Graphically it represents the area below the function, as picture shows below.

In many numerical methods like FEM, FVM, etc it is necessary to calculate integral numerically. There are several methods but the Gauss Quadrature is the most used one and popular.

Gaus Quadrature method of integration is based on the fact that if we make transformation of the function between interval in to another function on interval we can calculate approximate value of the integral on very simple way.

The picture below shows how we transform out starting function in to function defined on interval .

If we make such a transformation numerically we can calculate integral on the following way:

The summation function is called the Legendre-Gauss quadrature rule because the abscissa in the Gauss quadrature function for are defined as the roots of the Legendre polynomial for :

,

with the weights coming from the following function:

.

Solution for and  are based on number of points which can be calculated by using several method implementation you can find on internet. In this demo values for  and    for n=2 to 64  is taken from this link.

C# Implementation

Gauss Quadrature 1D numeric integrator demo is shown on the picture below.

It is a  small Windows Forms application by which you can calculate almost any definite integral. Demo uses Flee – Fast Lightweight Expression Evaluator for evaluating string in to math expression. As picture shows write down any function with x as an independent variable, define interval of integration, specify number of Gaussian points and press Calculate button.

Demo contains external file for Gaussian points and weights for n=2 to 64 points. So you can integrate any function with up to 64 Gaussian points which represent very large range of numerical application.

1. The first step in implementation was to implement loading Gaussian points from external files. The following implementation shows loading points to .NET object.

private List LoadGaus1DValues()
{
    gaus1DValues = new List();
    gaus1DValues.Add(new Helpers());
    string buffer = "";
    // open  file and retrieve the content
    using (StreamReader reader = System.IO.File.OpenText("GaussValues_n1_64_1d.csv"))
    {
        //read data in to buffer
        buffer = reader.ReadToEnd();
        reader.DiscardBufferedData();
        //reader.Close();
    }

    //define the lines from file
    var lines = (from l in buffer.Split(new char[] { '\n' }, StringSplitOptions.RemoveEmptyEntries)
                    where l[0] != '!' && l[0] != '\r'
                    select l.IndexOf('!') == -1 ? l : l.Remove(l.IndexOf('!'))
                    ).ToArray();

    int counter=1;
    for (int i=0; i<lines.Length; i++)
    {
        string line = lines[i].Replace("\r", "");
        string str = "n=" + counter.ToString();
        if (line == str)
        {
            i++;
            Helpers h = new Helpers();
            h.n = counter;
            h.wi= new double[counter];
            h.xi = new double[counter];
            for (int j = 0; j < counter; j++)
            {
                var lns= lines[i].Split(new char[] { '\t' });
                h.wi[j] = double.Parse(lns[1], CultureInfo.InvariantCulture);
                var s = lns[2].Replace("\r","");
                h.xi[j] = double.Parse(s, CultureInfo.InvariantCulture);
                if(j+1<counter)
                    i++;
            }
            counter++;
            gaus1DValues.Add(h);
        }
    }
    return gaus1DValues;
}

As can be seen from the code above, all values are stored in list of object of  Helper class which contains Gaussian values for each integration point. When the user click on Calculate button, the following code is going to be executed:

try
{
    var funstr = textFunction.Text;
	function = context.CompileGeneric(funstr);
}
catch (Exception ex)
{
    MessageBox.Show("Error in fun definition. ERROR: "+ex.Message);
    label5.Text = "I=n/a";
	return;
}

//perform integration
double retVal= Integrate1D(a,b,n);
label5.Text = string.Format("I={0:0.0000000000}", retVal);

Integrate1D(a,b,n); – method is implemented on the following way:

private double Integrate1D(double a, double b, int n)
{
    double I = 0;
    for (int i = 0; i < n; i++)
    {
        var ksi = 0.5 * (a + b) + 0.5 * (b - a) * gaus1DValues[n].xi[i];//f(1/2(a+b)+1/2*(a-b)xi)
        var dx = 0.5 * (b - a);//det|J|
        register.x = ksi;
        var f = function.Evaluate();
        //intgreation
        I += f * dx * gaus1DValues[n].wi[i];
    }

    return I;
}

As you can see this is very easy implementation,by using great expression evaluator and predefined Gaussian points.

The next blog post will be implementation of integral for 2D problems similar we can see on classic FEA engineering problems of stress, fracture mechanic etc..