Monthly Archives: July 2013

Gauss Numeric Integrator – my new open source project for numerical integration


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

C# implementation of numerical integration by Gaussian Quadrature of 1D function


Update: Gauss Numeric Integrator as open source project hosted on codeplex:

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

As you know definite integral can be expressed as:

\int_{a}^{b}{f(x)} dx.

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 f(x) between interval (a,b) in to another function \phi ( \xi ) on interval (-1,1) 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 (-1,1).

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

\int_{a}^{b}{f(x)}dx = \frac{a-b}{2} \int_{-1}^{1}{f(\frac{a-b}{2} x_i+ \frac{a+b}{2})}dx = \frac{a-b}{2}\sum_{i=1}^{n} w_if(\frac{a-b}{2} x_i+ \frac{a+b}{2})

The summation function is called the Legendre-Gauss quadrature rule because the abscissa x_i in the Gauss quadrature function for [-1,1] are defined as the roots of the Legendre polynomial for n:

P_n(x)= \frac{1}{2 \pi i} \oint (1-2tx+r^2)^{-1/2}t^{-n-1},

with the weights w_i coming from the following function:

w_i= \frac{-2} {(1-x_i^2)[P'_n(x_i))]^2} .

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

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

GPdotNET v3.0 is out


I am very proud to announce GPdotNET v3.0, the new version which brings 3 new solvers based on Genetic Algorithm and several new features. In previous posts I was  writing about those new features:

1. TSP and vector based chromosome implementation in GPdotNET

2. Exporting GP Model in to Wolfram Mathematica

3. Matrix based chromosome implementation

4. Introduction to Assignment and transportation problems based on Genetic Algorithm

The source code and Click Once installation are uploaded in codeplex site, so go and grab the new version. hppt://gpdotnet.codeplex.com

I am asking for feedback and bug reports. If you find something strange post a comment here on anywhere in blog site or codeplex site. I will try to answer as quickly as possible.

Here are screenshots with new Assignment and Transportation solvers in GPdotNET:


Loaded  training data of Assignment problem

Simulation Tab page with optimal result of Assignment problem

Loaded  training data of Transportation problem


Simulation Tab page with optimal result of Transportation problem

Get every new post delivered to your Inbox.

Join 672 other followers