Create lgwt.m

Author: Numerically-Stable

Numerical Methods/Gauss Legendre Quadrature/lgwt.m

@@ -0,0 +1,36 @@
+function [x, w] = lgwt(N, a, b)
+    % lgwt.m - Compute Gauss Legendre nodes and weights
+    % From 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
+    Lp = zeros(N1,1);
+
+    % Iterate using Newton's method
+    y0 = 2;
+    while max(abs(y - y0)) > eps
+        L(:,1) = 1;
+        L(:,2) = y;
+
+        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;
+    w = (b - a)./((1 - y.^2).*Lp.^2)*(N2/N1)^2;
+end