Gauss Quadrature in C#

[Tsotne]

Hello again! I need a help in identifying the problem why I don't get the correct answer in the following algorithm. I have defined the Gauss Quadrature algorithm and I always get the integral between [0,1] interval whatever the function I use.

static double Iteration(double[] qsi, double[] w, double delta){
                double[] qsi_5 = { 0, 0.5384693101056830910363144, 0.9061798459386639927976269 };
                double[] w_5 = { 0.5688888888888888888888889,
                                   0.4786286704993664680412915,
                                   0.2369268850561890875142640 };

                return Iteration(qsi_5, w_5, delta);
}

static double Iteration(double[] qsi, double[] w, double delta)
        {
            double S = 0;
            int n = qsi.Length;

            if (qsi[0] == 0)
            {
                for (int i = 1; i < n; i++)
                {
                    S += w[i] * F(delta * qsi[i] + delta);
                    S += w[i] * F(-delta * qsi[i] + delta);
                }
                S += w[0] * F(delta);
                S *= delta;
            }
            else
            {
                for (int i = 0; i < n; i++)
                {
                    S += w[i] * F(delta * qsi[i] + delta);
                    S += w[i] * F(-delta * qsi[i] + delta);
                }
                S *= delta;
            }
            return S;
        }

Here assume that the F function is not delegated but simply defined and invoked in the method.

Thank you

 

[cgorac]

The procedure above will produce correct answer, but for the fixed integration interval [0,2*delta]. I guess you may wish to change the code, in order to support integration over the arbitrary interval [a,b], to something as follows:

class Test
{
        static double F(double x)
        {
                return System.Math.Sin(x);
        }

        static double GaussianQuadrature(double[] qsi, double[] w, double a, double b)
        {
                double delta = (b - a) / 2;
                double mid = (a + b) / 2;
                double S = 0;
                int n = qsi.Length;

                if (qsi[0] == 0) {
                        for (int i = 1; i < n; i++) {
                                S += w[i] * F(delta * qsi[i] + mid);
                                S += w[i] * F(-delta * qsi[i] + mid);
                        }
                        S += w[0] * F(mid);
                        S *= delta;
                }
                else{
                        for (int i = 0; i < n; i++)  {
                                S += w[i] * F(delta * qsi[i] + mid);
                                S += w[i] * F(-delta * qsi[i] + mid);
                        }
                        S *= delta;
                }
                return S;
        }

        static void Main()
        {
                double[] qsi_5 = {
                        0,
                        0.5384693101056830910363144,
                        0.9061798459386639927976269
                };
                double[] w_5 = {
                        0.5688888888888888888888889,
                        0.4786286704993664680412915,
                        0.2369268850561890875142640
                };
                System.Console.WriteLine(GaussianQuadrature(qsi_5, w_5, System.Math.PI / 6, System.Math.PI / 2));
        }
}

 

[Tsotne]

Thank you very much @cgorac for your response. I've been waiting for a response so long I even forgot this thread. Thanks again.

I think I didn't exhaustively explained the above code I wrote. The qsi_5 and w_5 actually should contain 5 members but since the members except 0 are symmetric to each other in regard to sign, I just included the minus sign and added again in the iteration code above.

S += w[i] * F(-delta * qsi[i] + mid);

I think you guessed this since you responded. Anyway it was my fault not to explain it fully.
Are you sure the code you generated produces the correct answer???

Thank you again

 

[cgorac]

Well I haven't run any kind of battery of tests over this code to be sure (you should certainly do this in case you intend to deploy alike code in any kind of production environment: use either some well-known Gaussian quadrature code, or a tool like Mathematica, to calculate integrals of number of functions over number of intervals, and compare results with the results produced by your code), but - yes, I reasonably believe that with the change indicated the code should calculate Gaussian quadratures correctly. But the real question here is: do you understand Gaussian quadrature fully? If so, it should not be hard to you to see the problem with your initial code, and how the change above should fix it.

 

[Tsotne]

Thanks.

do you understand Gaussian quadrature fully? If so, it should not be hard to you to see the problem with your initial code, and how the change above should fix it.

Yes I understood what change you made and what you meant in it. But exactly as you say, I have defined some other algorithms like Gauss_Kronod, Simpson, etc. and when I compare it I still get different answers. I think the problem is not with the algorithm formula or programming constructs themselves. I'll fix it later when I get home. Thanks again for your help

 

[Tsotne]

I have also defined the Simpson algorithm to evaluate integrals. It has a good precision and is quite fast with small intervals as it is considered to be while using this algorithms. But the problem is that it doesn't at all evaluate the integrals of the functions like this:

Don't know why. And when I put the integration bounds as -4.999999999999 and 4.999999999999 it gives the approximate answer. Why cannot it evaluate the integrand over [-5,5]? Thank you very much. Here is the code:

static double Simpson(double a, double b)
        {
            double A = 0;
            double A0 = 0;
            double A1 = 0;
            double sum = 0;
            int n = 100;

            do
            {
                sum = 0;
                double deltaX = (b - a) / n;
                double[] Y = new double[n + 1];

                for (int i = 0; i <= n; i++)
                {
                    Y[i] = F(a + i * deltaX);
                }

                for (int i = 1; i <= n - 1; i += 2)
                {
                    sum += 4 * Y[i];
                }

                for (int j = 2; j <= n - 2; j += 2)
                {
                    sum += 2 * Y[j];
                }

                A = (deltaX / 3) * (Y[0] + sum + Y[n]);
                n *= 2;
                A0 = A1;
                A1 = A;
            } while (Math.Abs(A1 - A0) > 10e-10);

            return A;
        }

For some purposes consider F(x) is a static function defined somewhere else.