In the book of "Implementing Spectral Methods for Partial Differential Equations" by David A. Kopriva it is written on page no. 34 (ch. # 1): For problems with Legendre weighted integrals, the abscissas and weights for the Gauss-Lobatto rule are calculated as: The abscissas (nodes) are
$x_j$ = +1, -1, Zeros of $L'_N(x)$
the weights are
$w_j = \frac{2}{N(N+1)}\frac{1}{[L_N(x_j)]^2}$ where j = 0, 1,...,N
by using these formulas kindly help me calculating nodes and weights for N = 16? how to calculate Zeros of $L'_N(x)$?
The Legendre Polynomials have a recursive formula: $(n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)$ (Here $P_n(x) = L_n(x)$ in your notation above). You can also use an explicit formula to find them: $P_{n}(x) ={\frac {1}{2^{n}}}\sum _{k=0}^{n}{n \choose k}^{2}(x-1)^{n-k}(x+1)^{k}$. Mathematica will find the polynomial with the command LegendreP[n,x]. The polynomial you are looking for is $$P_{16}(x) = (1/32768)(6435 - 875160 x^2 + 19399380 x^4 - 162954792 x^6 + 669278610 x^8 - 1487285800 x^{10} + 1825305300 x^{12} - 1163381400 x^{14} + 300540195 x^{16}).$$ The derivative is $$P'_{16}(x) = (1/32768)(-1750320 x + 77597520 x^3 - 977728752 x^5 + 5354228880 x^7 - 14872858000 x^9 + 21903663600 x^{11} - 16287339600 x^{13} + 4808643120 x^{15})$$ which is a 15th order polynomial with 15 real roots. I have no hope of find the roots of this polynomial by hand, but a computer can find approximations (with Mathematica, try Solve[D[LegendreP[16,x],x] == 0.0]).
You'll get 15 roots: $[-0.973132, -0.91088, -0.815696, -0.691029, -0.541385, -0.372174, -0.189512, 0, 0.189512, 0.372174, 0.541385, 0.691029, 0.815696, 0.91088, 0.973132]$
Add $-1$ to the beginning of this list, and $1$ to the end of it, and you have your nodes $x_i$. Using $P_{16}(x)$, you can now find your weights. For $x_0 = -1$ and $x_{16} = 1$, you'll use the weights $$w_0 = w_{16} = \frac{2}{n(n+1)} = \frac{2}{(16)(17)}.$$ For $x_1 = -0.973132$, for example, you'll use the weight $$w_1 = \frac{2}{n(n+1)P_n(x_1)} = \frac{2}{(16)(17)P_n(-0.973132)} = -0.0181744.$$
Then, to do the Gauss-Lobatto integration rule for $\int_{-1}^{1} f(x)dx$ you compute $\sum_{i=0}^n w_i f(x_i)$.
The casio website computes the nodes and weights for you, but their formula starts at $n=1$ instead of $n=0$, so it is off by one of yours, and you'll want to compute for $n=17$ instead on that website.