# Gauss Legendre quadrature problem with Legendre polynomials composed with square root

Let $$P_n$$ be the orthogonal Legendre polynomial with a degree of $$n$$, meaning it satisfies the following recursive formula: $$(n+1)P_{n+1}(x)-(2n+1)xP_n(x)+nP_{n-1}(x)=0$$ where $$P_0(x) = 1$$ and $$P_1(x) = x$$. Let $$w(x)=1$$ be the weight function, and the interval $$[-1,1]$$. $$I_n(f)=\sum_{k=1}^nw_k f(x_k)$$ is the integration formula. For $$n\geq1$$ we define $$q_n(y) =\frac{P_{2n+1}(\sqrt{y})}{\sqrt{y}}.$$

Problem: prove that $$q_n$$ is a polynomial of degree $$n$$, determine the orthogonality relations, weight function (for the dot product $$\int_{-1}^1w(x)f(x)g(x)dx$$) and the interval where they form a set of orthogonal polynomials. Express the weights and nodes of the formula $$I_n^q =\sum_{k=1}^n v_k f(y_k)$$ using $$w_k$$ and $$x_k$$ from the integration formula $$I_{2n+1}$$.

My attempt:

I proved that $$q_n$$ is a polynomial of degree $$n$$ using induction and the fact that Legendre polynomial of degree $$n$$ is even/odd if $$n$$ is even/odd.

Next, I determined that they are orthogonal with the weight function being $$x \mapsto \sqrt{x}$$ and the interval $$[0,1]$$, by using the fact that Legendre polynomials are orthogonal and some substitutions.

Now, given that the nodes in the formula $$I_n$$ are roots of the n-th orthogonal polynomial we're using, when we look at the formula $$I_{2n+1}$$ we're looking for roots of $$P_{2n+1}$$. Seeing as it's a polynomial of degree $$2n+1$$ that is odd, it must have $$n$$ positive roots (let's denote them by $$x_1,\ldots x_n$$, $$n$$ negative roots and $$0$$. (there's also a theorem that states that the roots of these polynomials are all different)

Now, I think I should take $$y_k = x_k^2$$, $$k = 1,\ldots,n$$.

But I'm not sure how to express $$v_k$$ using $$w_k$$.

Any hints would be appreciated!

I don't see an elegant way to do this so we will apply brute force starting from this formula $$w_k=\frac{a_n}{a_{n-1}}\frac{\int_a^bw(x)p_{n-1 }(x)^2dx}{p_n^{\prime}(x_k)p_{n-1}(x_k)}$$ In our case we write it as $$v_k=\frac{b_m}{b_{m-1}}\frac{\int_0^1\sqrt y\,q_{m-1}(y)^2dy}{q_m^{\prime}(y_k)q_{m-1}(y_k)}$$ Now $$b_m$$ is the leading coefficient of $$q_m(y)=\frac{P_{2m+1}(\sqrt y)}{\sqrt y}$$ so $$b_m=a_{2m+1}$$, the leading coefficient of $$P_{2m+1}(x)$$. This is a good start! Now, $$b_{m-1}=a_{2m-1}$$: not quite what we want but hopefully it will work itself out. We have $$q_m^{\prime}(y_k)=\frac{P_{2m+1}^{\prime}(\sqrt{y_k})}{\sqrt{y_k}}\frac1{2\sqrt{y_k}}-\frac12\frac{P_{2m+1}(y_k)}{y_k^{3/2}}=\frac{P_{2m+1}^{\prime}(x_k)}{2x_k^2}$$ $$q_{m-1}(y_k)=\frac{P_{2m-1}(\sqrt{y_k})}{\sqrt{y_k}}=\frac{P_{2m-1}(x_k)}{x_k}$$ We can adjust this equation by taking $$n=2m$$ and $$x=x_k$$ the the three term recurrence relation to get $$(2m+1)P_{2m+1}(x_k)-(4m+1)x_kP_{2m}(x_k)+2mP_{2m-1}(x_k)=0$$ Solve for $$P_{2m-1}(x_k)$$ and apply to the previous equation to get $$q_{m-1}(y_k)=\frac{4m+1}{2m}P_{2m}(x_k)$$ Which is the form we want. Now, \begin{align}\int_0^1\sqrt y\,q_{m-1}(y)^2dy&=\int_0^1\sqrt y\,\frac{P_{2m-1}(\sqrt y)^2}{y}dy=\int_0^1x\frac{P_{2m-1}(x)^2}{x^2}\cdot2x\,dx\\ &=2\int_0^1P_{2m-1}(x)^2dx=2\int_0^{-1}P_{2m-1}(-u)^2(-du)\\ &=2\int_{-1}^0\left(-P_{2m-1}(u)\right)^2du=\int_{-1}^1P_{2m-1}(x)^2dx\end{align} We go bck to the three term recurrence relation with $$n=2m$$, multiply by $$P_{2m-1}(x)$$ and integerate from $$-1$$ to $$1$$ with respect to $$x$$ to get $$(2m+1)\int_{-1}^1P_{2m+1}(x)P_{2m-1}(x)dx-(4m+1)\int_{-1}^1P_{2m}(x)\left[xP_{2m-1}(x)\right]dx+2m\int_{-1}^1P_{2m-1}(x)^2dx=0$$ The first integral above is $$0$$ by the orthogonality of the Legendre polynomials and $$xP_{2m-1}(x)=\frac{a_{2m-1}}{a_{2m}}P_{2m}(x)+R_{2m-1}(x)$$ By the division algorithm where $$R_{2m-1}(x)$$ is a polynomial in $$x$$ of degree at most $$2m-1$$. Thus by the orthogonality of the Legendre polynomials, \begin{align}\int_{-1}^1P_{2m}(x)\left[xP_{2m-1}(x)\right]dx&=\int_{-1}^1P_{2m}(x)\left[\frac{a_{2m-1}}{a_{2m}}P_{2m}(x)+R_{2m-1}(x)\right]dx\\ &=\frac{a_{2m-1}}{a_{2m}}\int_{-1}^1P_{2m}(x)^2dx\end{align} So we have arrived at $$-(4m+1)\frac{a_{2m-1}}{a_{2m}}\int_{-1}^1P_{2m}(x)^2dx+2m\int_{-1}^1P_{2m-1}(x)^2dx=0$$ And we can use this result to achieve $$\int_0^1\sqrt y\,q_{m-1}(y)^2dy=\frac{(4m+1)a_{2m-1}}{2ma_{2m}}\int_{-1}^1P_{2m}(x)^2dx$$ So now we have \begin{align}v_k&=\frac{a_{2m+1}}{a_{2m-1}}\frac{\frac{(4m+1)a_{2m-1}}{2ma_{2m}}\int_{-1}^1P_{2m}(x)^2dx}{\frac{P_{2m+1}^{\prime}(x_k)}{2x_k^2}\frac{(4m+1)P_{2m}(x_k)}{2m}}\\ &=2x_k^2\frac{a_{2m+1}}{a_{2m}}\frac{\int_{-1}^1P_{2m+1}(x)^2dx}{P_{2m+1}^{\prime}(x_k)P_{2m}(x_k)}\\ &=2x_k^2w_k\end{align} Mmm... now that we have an answer, it looks like an elegant proof should be possible, but I don't see it just now. It's getting late so I leave the reader with the above brute force approach.
EDIT: It's late but not that late. We start with the most horrible possible way of computing the weights $$v_k$$: $$I_j=\int_0^1\sqrt y\,y^jdy=\sum_{k=1}^mv_ky_k^j=\sum_{k=1}^mv_kx_k^{2j}$$ For $$0\le j\le m-1$$, where we have applied our Gauss-Jacobi quadrature formula, exact for polynomials of degree at most $$2m-1$$, to find $$m$$ linear equations for its $$m$$ weights. Since the matrix of coefficients is a Vandermonde matrix, the weights are uniquely determined by this system.
Now we transform variables first to get \begin{align}I_j&=\int_0^1\sqrt y\,y^jdy=\int_0^1x\cdot x^{2j}\cdot 2x\,dx=\int_{-1}^1x^{2j+2}dx\\ &=2\sum_{k=1}^mw_kx_k^{2j+2}=\sum_{k=1}^m(2x_k^2w_k)x_k^{2j}\end{align} For $$0\le j\le m-1$$. We showed in our complicated proof how to fold the integral over the origin to change the interval to $$[-1,1]$$ and the Gauss-Legendre quadrature formula was again exact for polynomials of degree at most $$4m+1$$. We know already that the zeros of the Legendre polynomial of degree $$2m+1$$ are symmetrically disposed about the origin and that the weight corresponding to $$-x_k$$ is the same as that corresponding to $$x_k$$, due to the symmetry of the weight function and interval in Gauss-Legendre quadrature about the origin. Also the weight for the zero node doesn't enter in because it always gets multiplied by zero in every equation.
Thus we see that $$v_k$$ and $$2x_k^2w_k$$ satisfy the same nonsingular system of $$m$$ equations in $$m$$ unknowns, so they are the same.