How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the D programming language

Published on 12 May 2024 09:40 PM
#D

How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the D programming language

Table of Contents

Problem Statement

For this, we first need to calculate the nodes and the weights, but after we have them, we can reuse them for numerious integral evaluations, which greatly speeds up the calculation compared to more simple numerical integration methods.

Task description Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function

f ( x )

{\displaystyle f(x)}

, but by applying an n-point Gauss-Legendre quadrature rule, as described here, for example. The input values should be an function f to integrate, the bounds of the integration interval a and b, and the number of gaussian evaluation points n. An reference implementation in Common Lisp is provided for comparison. To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the D programming language

Source code in the d programming language

import std.stdio, std.math;

immutable struct GaussLegendreQuadrature(size_t N, FP=double,
                                         size_t NBITS=50) {
    immutable static double[N] lroots, weight;
    alias FP[N + 1][N + 1] CoefMat;

    pure nothrow @safe @nogc static this() {
        static FP legendreEval(in ref FP[N + 1][N + 1] lcoef,
                               in int n, in FP x) pure nothrow {
            FP s = lcoef[n][n];
            foreach_reverse (immutable i; 1 .. n+1)
                s = s * x + lcoef[n][i - 1];
            return s;
        }

        static FP legendreDiff(in ref CoefMat lcoef,
                               in int n, in FP x)
        pure nothrow @safe @nogc {
            return n * (x * legendreEval(lcoef, n, x) -
                        legendreEval(lcoef, n - 1, x)) /
                   (x ^^ 2 - 1);
        }

        CoefMat lcoef = 0.0;
        legendreCoefInit(/*ref*/ lcoef);

        // Legendre roots:
        foreach (immutable i; 1 .. N + 1) {
            FP x = cos(PI * (i - 0.25) / (N + 0.5));
            FP x1;
            do {
                x1 = x;
                x -= legendreEval(lcoef, N, x) /
                     legendreDiff(lcoef, N, x);
            } while (feqrel(x, x1) < NBITS);
            lroots[i - 1] = x;
            x1 = legendreDiff(lcoef, N, x);
            weight[i - 1] = 2 / ((1 - x ^^ 2) * (x1 ^^ 2));
        }
    }

    static private void legendreCoefInit(ref CoefMat lcoef)
    pure nothrow @safe @nogc {
        lcoef[0][0] = lcoef[1][1] = 1;
        foreach (immutable int n; 2 .. N + 1) { // n must be signed.
            lcoef[n][0] = -(n - 1) * lcoef[n - 2][0] / n;
            foreach (immutable i; 1 .. n + 1)
                lcoef[n][i] = ((2 * n - 1) * lcoef[n - 1][i - 1] -
                               (n - 1) * lcoef[n - 2][i]) / n;
        }
    }

    static public FP integrate(in FP function(/*in*/ FP x) pure nothrow @safe @nogc f,
                               in FP a, in FP b)
    pure nothrow @safe @nogc {
        immutable FP c1 = (b - a) / 2;
        immutable FP c2 = (b + a) / 2;
        FP sum = 0.0;
        foreach (immutable i; 0 .. N)
            sum += weight[i] * f(c1 * lroots[i] + c2);
        return c1 * sum;
    }
}

void main() {
    GaussLegendreQuadrature!(5, real) glq;
    writeln("Roots:  ", glq.lroots);
    writeln("Weight: ", glq.weight);
    writefln("Integrating exp(x) over [-3, 3]: %10.12f",
             glq.integrate(&exp, -3, 3));
    writefln("Compred to actual:               %10.12f",
             3.0.exp - exp(-3.0));
}

You may also check:How to resolve the algorithm User input/Graphical step by step in the Ruby programming language
You may also check:How to resolve the algorithm Roman numerals/Encode step by step in the Red programming language
You may also check:How to resolve the algorithm Long primes step by step in the Quackery programming language
You may also check:How to resolve the algorithm Count occurrences of a substring step by step in the Tcl programming language
You may also check:How to resolve the algorithm Cistercian numerals step by step in the C++ programming language