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


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

transformation1

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

transformation

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

C# Implementation

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

transformation2

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

About Bahrudin Hrnjica

PhD in Mechanical Engineering, Microsoft MVP for .NET. Likes .NET, Math, Mechanical Engineering, Evolutionary Algorithms, Blogging.

Posted on 12/07/2013, in .NET, C#, Numeric Calculation and tagged , , , , . Bookmark the permalink. 3 Comments.

  1. I read a lot of interesting articles here. Probably you spend
    a lot of time writing, i know how to save you a lot of
    work, there is an online tool that creates readable, google
    friendly posts in minutes, just search in google – k2seotips
    unlimited content

  2. Your description of this issue is very clear, so that we can understand,.You are very excellent, thank you for your article ahout c#implementation,it’s very useful to me.I am going to read your next blog about 2D problems.

  1. Pingback: Gauss Numeric Integrator – my new open source project for numerical integration | Bahrudin Hrnjica Blog

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s