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

Published on 12 May 2024 09:40 PM
#Pl/i

How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the PL/I 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 PL/I programming language

Source code in the pl/i programming language

(subscriptrange, size, fofl):
Integration_Gauss: procedure options (main);

  declare (n, k) fixed binary;
  declare r(*,*) float (18) controlled;
  declare (z, a, b, exact) float (18);

  do n = 1 to 20;
    a = -3; b = 3;
    if allocation(r) > 0 then free r;
    allocate r(2, n); r = 0;
    call gaussquad(n, r);
    z = (b-a)/2 * sum(r(2,*) * exp((a+b)/2+r(1,*)*(b-a)/2));
    exact = exp(3.0q0)-exp(-3.0q0);
    put skip edit (n, z, z-exact) (f(5), f(25,16), e(15,2));
  end;

gaussquad: procedure(n, r);
/*declare n fixed binary, r(2, n) float (18);*/
  declare n fixed binary, r(2, *) float (18);/* corrected */
  declare pi float (18) value (4*atan(1.0q0));
  declare (x, f, df, dx) float (18);
  declare (i, iter, L) fixed binary;
  declare (p0(*), p1(*), tmp(*), tmp2(*)) float (18) controlled;
  
  allocate p0(1) initial (1);
  allocate p1(2) initial (1, 0);
  
  do k = 2 to n;
     allocate tmp(hbound(p1)+1); do L = 1 to hbound(p1); tmp(L) = p1(L); end; tmp(L) = 0;
     allocate tmp2(hbound(p0)+2); tmp2(1), tmp2(2) = 0;
     do L = 1 to hbound(p0); tmp2(L+2) = p0(L); end;
     tmp = ((2*k-1)*tmp - (k-1)*tmp2)/k;
     free p0; allocate p0(hbound(p1)); p0 = p1;
     free p1; allocate p1(hbound(tmp)); p1 = tmp;
     free tmp, tmp2;
  end;
  do i = 1 to n;
    x = cos(pi*(i-0.25q0)/(n+0.5q0));
    do iter = 1 to 10;
      f = p1(1); df = 0;
      do k = 2 to hbound(p1);
        df = f + x*df;
        f  = p1(k) + x * f;
      end;
      dx =  f / df;
      x = x - dx;
      if abs(dx) < 10*epsilon(dx) then leave;
    end;
    r(1,i) = x;
    r(2,i) = 2/((1-x**2)*df**2);
  end;
  end gaussquad;
end Integration_Gauss;

You may also check:How to resolve the algorithm Loops/Foreach step by step in the AmigaE programming language
You may also check:How to resolve the algorithm Soundex step by step in the Java programming language
You may also check:How to resolve the algorithm Repeat a string step by step in the DWScript programming language
You may also check:How to resolve the algorithm Terminal control/Preserve screen step by step in the Scala programming language
You may also check:How to resolve the algorithm IBAN step by step in the UNIX Shell programming language