A Gaussian quadrature is used to approximate the following integral: $$ \int_{-1}^{1} f(x) dx \approx \sum_{i=1}^n w_i f(x_i). $$ Numerically I've found an interesting property of $x_i$ and $w_i$: if we split $[-1,1]$ interval into subintervals $I_i = [y_i, y_{i+1}]$ each of length $w_i$ (order matters) then each node $x_i$ belongs to corresponding interval $I_i$.
Formally speaking if we define $$ y_i = \sum_{k=1}^{i-1} w_k - 1 $$ then sequences $\{x_i\}_{i=1}^n$ and $\{y_i\}_{i=1}^{n+1}$ interleave.
Here's an example for $n=15$ (nodes are blue and the intervals endpoints' are marked with yellow)
Can this be proven for all $n$? I thought this could be proved using interleaving property of orthogonal polynomial roots, but $y_i$ do not seem to relate to any of those.
Edit Moreover, I've found the similar property holds for Gauss-Chebyshev quadrature. If we rewrite the quadrature formula as $$ \int_{-1}^1 f(x) \omega(x) dx \approx \sum_{i=1}^n f(x_i) \int_{y_i}^{y_{i+1}} \omega(x) dx, $$ and the $y_i$ are now obtained from solving $$ \int_{y_i}^{y_{i+1}} \omega(x) dx = w_i. $$