Menghitung bobot dan node pada Gauss Quadrature

Untuk menghitung integral berhingga dengan Gauss-Quadrature sebenarnya agak-agak rumit. Karena kita harus menempuh beberapa metode yang agak ruwet dan tidak mudah. Harus menyelesaikan sistem persamaan Linear dan lain-lain manipulasi aljabar yang membingungkan. Karena kita harus menghitung bobot dan node yang tergantung berapa nodenya. Tapi saya menemukan suatu source code di internet yang sudah menyelesaikan persoalan tersebut. Yakni

function [x,w]=lgwt(N,a,b)

% lgwt.m
%
% This script is for computing definite integrals using Legendre-Gauss 
% Quadrature. Computes the Legendre-Gauss nodes and weights  on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel - 02/25/2004
N=N-1;
N1=N+1; N2=N+2;

xu=linspace(-1,1,N1)';

% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);

% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);

% Derivative of LGVM
Lp=zeros(N1,N2);

% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method

y0=2;

% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
    
    
    L(:,1)=1;
    Lp(:,1)=0;
    
    L(:,2)=y;
    Lp(:,2)=1;
    
    for k=2:N1
        L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
    end
 
    Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);   
    
    y0=y;
    y=y0-L(:,N2)./Lp;
    
end

% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;      

% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;

Di mana versi python-nya adalah:

import numpy as np 

np.set_printoptions(precision=4)

# fungsi untuk menghitung node-node pada gauss-quadrature
def lgwt(N, a , b):
    N = N -1 
    N1 , N2= N + 1, N+2 
    xu = np.linspace(-1,1, N1) 
    y =  np.cos(( 2* np.transpose(np.arange(0,N+1)) + 1  )\
    * np.pi / (2 * N + 2))  + (0.27 / N1 ) * np.sin(np.pi * xu * N / N2) ; 
    
    L = np.zeros([N1, N2])
    y0 = 2 
    
    while np.max(np.abs(y - y0 )) > spacing(1):
        L[:,0] = 1        
        L[:,1] = y
        
        for k in range(1, N1):
            L[:,k+1] = ((2*(k+1)-1)*y *L[:,k]-(k)*L[:,k-1])/(k+1)
        
        Lp=(N2)*(L[:,N1-1]-y * L[:,N2-1])/(1-np.power(y,2))   

        y0 = y 
        y = y0 - L[:, N2-1] / Lp
    
    N2 = float(N2) 
    N1 = float(N1)
    x = (a*(1-y ) + b * (1+ y) ) / 2
    w = (b-a) / ((1 -  np.power(y ,2 ) ) * np.power(Lp, 2) ) \
        * np.power(N2/N1 , 2) 
    return (x, w) 
    
(a , b ) = lgwt(12,-5,5) # untuk ngetes
print  a 
print  b 

Potongan eksekusinya adalah

Anda bisa mengkonfirmasi di buku nilai-nilai di atas.

Kebetulan dasar teorinya agak ribet, tapi saya menemukan salah satunya di internet juga yakni sebuah tesis magister. Yang bisa anda download di:
https://drive.google.com/file/d/0B1irLqfPwjq0RFhHQVpvaUVuWlk/edit?usp=sharing