3

I am trying to understand the code given here, which calculates the Legendre-Gauss-Lobatto nodes via the Newton-Raphson method. These nodes are given by the zeros of $(1 - x^2)\,P'_{N}(x)$ where $P_N(x)$ is the Nth Legendre polynomial.

The code calculates these zeros via the Newton-Raphson method, which is used for finding the zeros of functions. For example, the zero of a function $f(x)$ may be found by iterating

$$x_{n+1} = x_{n} - \frac{f(x_n)}{f'(x_n)}$$

until the difference between $x_{n+1}$ and $x_n$ is sufficiently small. In the code, this update seems to be given by

$$x_{n+1} = x_n - \frac{x_n\,P_N(x_n) - P_{N-1}(x_n)}{N\,P_N(x_n)}.$$

The term on the right seems very similar to Bonnet's recursion formula

$$(x^2 - 1)\,P'_N(x) = N(x\,P_N(x) - P_{N-1}(x))$$

which mostly explains the numerator, but I am lost on the denominator.

If I differentiate Bonnet's recursion formula, I obtain

$$\frac{d}{dx} [(x^2 - 1)\,P'_N(x)] = N(x\,P'_N(x) + P_N(x) - P'_{N-1}(x))$$

which doesn't match the denominator. At the very least, I would expect there to be an $N^2$.

jltorchinsky
  • 163

1 Answers1

4

Okay, figured it out! The key is the differential equation definition of the Legendre Polynomials, which is

$$(x^2 - 1)\,P''_k(x) + 2\,x\,P'_k(x) - k\,(k+1)\,P_k(x) = 0$$

Differentiating $(x^2 - 1)\,P'_k(x)$ gives

$$\frac{d}{dx}\left[ (x^2 - 1)\,P'_k(x) \right] = (x^2 - 1)\,P''_k(x) + 2\,x\,P'_k(x) = k\,(k+1)\,P_k(x).$$

Hence, the second term on the right-hand side is

$$\frac{(x^2 - 1)\,P'_N(x)}{\frac{d}{dx}\left[ (x^2 - 1)\,P'_k(x) \right]} = \frac{x\,P_N(x) - P_{N-1}(x)}{(N + 1)\,P_N(x)}$$

with some indexing issues because Matlab indexes from 1 while the polynomials index from 0.

jltorchinsky
  • 163