Gauss Numeric Integrator Project

After short post about numerical integration in 1D, I have decided to start with open source project for numerical integration for any function in 1,2 or 3D.  If you are interesting in numeric analysis you can find a lot of stuff on the internet. For this blog post I prepared Demo example which extends the previous post by providing integration for 2D function. By using this demo you can integrate any 2D function on area defined by 3, 4 or 8 points. That means you can perform numerical integration for 3 kind of finite elements: triangle, quadrilateral 4-point or 8-point element. Beside element type you can choose up to 7 Gaussian integration points.

Implementation

The project contains one ClassLibrary and one Windows application (see picture below). Class library contains implementation related to numeric integration: Shape function, partial derivatives of shape function, Jacobina calculation for every element type as well as Integration process.

Windows Application implements Form for defining function, and other parameters. The project uses Flee – fast .Net expression evaluator for evaluating string in to math expression.

The Gauss Numeric Integrator can be downloaded from codeplex site at http://gauss.codeplex.com/

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..