We know,from generating function of Legendre polynomials,
$(1-2xz+z^2)^\frac{-1}{2}=\sum z^nP_n(x)$
Squaring both sides and writing the sum by separating the cross terms and "self-quadrature" terms,we get
$(1-2xz+z^2)^{-1}=\sum z^{2n}P_n^2(x)+2\sum z^{m+n}P_m(x)P_n(x)$
Integrating both sides between -1 and +1,we have,
$\int_{-1}^1\sum z^{2n}P_n^2(x)dx+\int_{-1}^1 2\sum z^{m+n}P_m(x)P_n(x)dx=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
Now,
$\int_{-1}^{1}P_n(x)P_m(x)dx=0$ due to orthogonality of Legendre polynomials and the 2nd integral on LHS of the above equation is wrt x,the terms containing m and n can be taken out,and therefore the 2nd integral vanishes.
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx+0=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
The integral on RHS,i.e.,
$
{\displaystyle\int}\dfrac{1}{-2zx+z^2+1}\,\mathrm{d}x $ can be solved as--
Substitute $u=-2zx+z^2+1$
$\dfrac{\mathrm{d}u}{\mathrm{d}x} = -2z$
$\Rightarrow \mathrm{d}x=-\dfrac{1}{2z}\,\mathrm{d}u$
The integral reduces to
$=-\class{steps-node}{\cssId{steps-node-1}{\dfrac{1}{2z}}}{\displaystyle\int}\dfrac{1}{u}\,\mathrm{d}u$
which is a standard integral.
Therefore the integral becomes
$=-\dfrac{\ln\left(u\right)}{2z}$
$=-\dfrac{\ln\left(-2zx+z^2+1\right)}{2z}$
$\Rightarrow \sum z^{2n}\int_{-1}^1P_n^2(x)dx=\frac{-1}{2z}[\ln(1-2xz+z^2)]_{-1}^1$
RHS=
$\frac{-1}{2z}\ln\frac{1-2z+z^2}{1+2z+z^2}=\frac{-1}{2z}\ln(\frac{1-z}{1+z})^2=\frac{1}{z}[\ln(1+z)-\ln(1-z)]$
Using Maclaurin Series of $\ln(1+z)$ and $\ln(1-z)$ and expanding,
RHS=
$\frac{1}{z}[(z-\frac{z^2}{2}+\frac{z^3}{3}-\frac{z^4}{4}+...)-(-z-\frac{z^2}{2}-\frac{z^3}{3}-\frac{z^4}{4}-...)]$
Observe that the terms containing even powers of x cancel and the odd terms gets doubled(due to presence of them in both the series indicated above).Taking the 1/z term inside the square brackets,
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx=2[1+\frac{z^2}{3}+\frac{z^4}{5}+..\frac{z^{2n}}{2n+1}.]$
Finally,equating the coefficients of $z^{2n}$ on both sides we have,
$\int_{-1}^{1} (P_n(x))^2 dx = \frac{2}{2n+1}$,which is the required relation.
What you’re doing is perhaps a little ‘unnatural’ given that, to me anyway, the more immediate approach, especially from a physics point of view is to consider the distance between two points in the plane and ask how we can get a simple series expansion i.e. let
$$\frac{1}{\sqrt{1-2xt+{{t}^{2}}}}=\sum\limits_{n=0}^{\infty }{{{P}_{n}}\left( x \right){{t}^{n}}}$$
and now tell me what ${{P}_{n}}\left( x \right)$ will make this happen. This I think is probably the origin of them. However starting with a second order ODE
$$\frac{d}{dz}\left( \left( 1-{{z}^{2}} \right)\frac{d}{dz}{{P}_{n}}\left( z \right) \right)+n\left( n+1 \right){{P}_{n}}\left( z \right)=0$$
and then ‘discovering’ the fact that they generate this inherently useful function, using your approach seems tricky. Someone out there might obviously see a very clean way of doing it but it’s not obvious to me anyway.
And one of the reasons I think this is might be that you’re utilising an approach where you’re generating a PDE
$$\frac{\partial }{\partial x}\left( \left( 1-{{x}^{2}} \right)\frac{\partial }{\partial x}g\left( x,t \right) \right)+\frac{\partial }{\partial t}{{t}^{2}}\frac{\partial }{\partial t}g\left( x,t \right)=0$$
Which you then need to solve with appropriate boundary conditions. Typically this is what sturm-liouville theory does because you take this, separate variables, and land at a linear ODE which have all these nice orthogonal properties. But in your case, this is exactly what you started with! So I think this is why you’re going in circles.
When this happens it’s telling you that you may need to take another angle. So in that regard, consider again
$$\frac{d}{dz}\left( \left( 1-{{z}^{2}} \right)\frac{d}{dz}{{y}_{n}}\left( z \right) \right)+n\left( n+1 \right){{y}_{n}}\left( z \right)=0$$
And lets sketch out a pathway to the generating function. Note $z=\pm 1$ are regular singular points and so the solution of this equation is analytic in the domain $\left| z \right|<1$. More so those solutions contain polynomials when n is an integer. To see this let ${{y}_{n}}\left( z \right)=\sum\limits_{m=0}^{\infty }{{{p}_{m}}{{z}^{m}}}$ and so
$$\sum\limits_{m=0}^{\infty }{m\left( m-1 \right){{p}_{m}}{{z}^{m-2}}}-2z\sum\limits_{m=0}^{\infty }{m{{p}_{m}}{{z}^{m-1}}}\\-\sum\limits_{m=0}^{\infty }{m\left( m-1 \right){{p}_{m}}{{z}^{m}}}+n\left( n+1 \right)\sum\limits_{m=0}^{\infty }{{{p}_{m}}{{z}^{m}}}=0$$
Or
$${{p}_{m+2}}=\frac{\left( m+1 \right)m-\left( n+1 \right)n}{\left( m+2 \right)\left( m+1 \right)}{{p}_{m}}$$
Next consider n to be a positive integer, and for now assume it to be even. We have initial conditions on the recurrence above being ${{p}_{0}}=a,\,\,{{p}_{1}}=b$ and so if we break the recurrence over even m and odd m, we obtain a solution in the form
$${{y}_{n}}\left( z \right)=\sum\limits_{m=even}^{{}}{{{p}_{m}}{{z}^{m}}}+\sum\limits_{m=odd}^{{}}{{{p}_{m}}{{z}^{m}}}$$
At some point, since n is even, we will get for even $m=n$ a coefficient which is zero and hence all subsequent even coefficients will also be zero, i.e. ${{p}_{even\,\,m>n}}=0$. So our solution becomes
$${{y}_{n}}\left( z \right)=\sum\limits_{m=0}^{n/2}{{{p}_{2m}}{{z}^{2m}}}+\sum\limits_{m=0}^{\infty }{{{p}_{2m+1}}{{z}^{2m+1}}}={{P}_{n}}\left( z \right)+H\left( z \right)$$
Where ${{P}_{n}}\left( z \right)$ is a polynomial of degree n, and H is an analytic function represented by an infinite series converging for $\left| z \right|<1$. A similar situation occurs when n is an odd positive integer. Now if we require our solution to be bounded at $\left| z \right|=1$ we effectively want to make $H\left( z \right)$ disappear (because you can show it diverges). This is easily achieved by choosing conditions which make the even or odd coefficients disappear (i.e. for even $n$ let $b=0$ and for odd $n$ let $a=0$). In this instance our solution to the DE is simply ${{y}_{n}}\left( z \right)={{P}_{n}}\left( z \right)$ which is assured to be analytic for $\left| z \right|\le 1$. These are the solutions we will focus on.
So far we have therefore shown that
$$\frac{d}{dz}\left( \left( 1-{{z}^{2}} \right)\frac{d}{dz}{{P}_{n}}\left( z \right) \right)+n\left( n+1 \right){{P}_{n}}\left( z \right)=0$$
where ${{P}_{n}}\left( z \right)$are n-th order polynomials. Now observe that from the above ODE we may write
$${{P}_{n}}\left( z \right)=-\frac{1}{n\left( n+1 \right)}\frac{d}{dz}\left( \left( 1-{{z}^{2}} \right)\frac{d}{dz}{{P}_{n}}\left( z \right) \right)$$
Which if we substitute for ${{P}_{n}}\left( z \right)$ again we have
$${{P}_{n}}\left( z \right)=\frac{1}{{{n}^{2}}{{\left( n+1 \right)}^{2}}}\frac{d}{dz}\left( \left( 1-{{z}^{2}} \right)\frac{{{d}^{2}}}{d{{z}^{2}}}\left( \left( 1-{{z}^{2}} \right)\frac{d}{dz}{{P}_{n}}\left( z \right) \right) \right)$$
And supposing we could carrying this on indefinitely we’d have a string of derivatives of the polynomial $\left( 1-{{z}^{2}} \right)$ . So let’s examine this and consider
$$f={{\left( 1-{{z}^{2}} \right)}^{n}}\Rightarrow f'+2nz{{\left( 1-{{z}^{2}} \right)}^{n-1}}=0\Rightarrow \left( 1-{{z}^{2}} \right)f'+2nzf=0$$
And we’re almost at our original DE. Now consider
$$\frac{{{d}^{k}}}{d{{z}^{k}}}\left( 1-{{z}^{2}} \right)f'=\sum\limits_{m=0}^{k}{\left( \begin{matrix}
k \\
m \\
\end{matrix} \right){{\left( 1-{{z}^{2}} \right)}^{\left( m \right)}}{{f}^{\left( k-m+1 \right)}}}$$
Where we’ve used Liebnitz formula
$$\frac{{{d}^{k}}}{d{{z}^{k}}}\left( wy' \right)=\sum\limits_{m=0}^{k}{\left( \begin{matrix}
k \\
m \\
\end{matrix} \right){{w}^{\left( m \right)}}{{y}^{\left( k-m+1 \right)}}}$$
Now since $w\left( z \right)=1-{{z}^{2}}$ is a second order polynomial, there can be only three terms and we therefore have
$$\frac{{{d}^{k}}}{d{{z}^{k}}}\left( 1-{{z}^{2}} \right)f'=\left( 1-{{z}^{2}} \right){{f}^{\left( k+1 \right)}}-2kz{{f}^{\left( k \right)}}-k\left( k-1 \right){{f}^{\left( k-1 \right)}}$$
Hence
$$\frac{{{d}^{k}}}{d{{z}^{k}}}\left\{ \left( 1-{{z}^{2}} \right)f'+2nzf \right\}=0\\\Rightarrow \left( 1-{{z}^{2}} \right){{f}^{\left( k+1 \right)}}-2kz{{f}^{\left( k \right)}}-k\left( k-1 \right){{f}^{\left( k-1 \right)}}+\frac{{{d}^{k}}}{d{{z}^{k}}}2nzf=0$$
Applying the rule again
$$\left( 1-{{z}^{2}} \right){{f}^{\left( k+1 \right)}}-2kz{{f}^{\left( k \right)}}-k\left( k-1 \right){{f}^{\left( k-1 \right)}}+2nz{{f}^{\left( k \right)}}+2nk{{f}^{\left( k-1 \right)}}=0$$
Let $k=n+1$ we have finally
$$\left( 1-{{z}^{2}} \right)y''-2zy'+n\left( n+1 \right)y=0$$
where $y={{f}^{\left( n \right)}}$. This is the DE we started with. So we now have that the polynomials might be able to be represented by the Rodrigues formula
$${{P}_{n}}\left( z \right)=A\frac{{{d}^{n}}}{d{{z}^{n}}}{{\left( 1-{{z}^{2}} \right)}^{n}}$$
for some constant A. The constant comes from the boundary condition we enforce upon say ${{P}_{n}}\left( 1 \right)$ which traditionally is ${{P}_{n}}\left( 1 \right)=1$. Doing so we find then
$${{P}_{n}}\left( z \right)=\frac{{{\left( -1 \right)}^{n}}}{{{2}^{n}}n!}\frac{{{d}^{n}}}{d{{z}^{n}}}{{\left( 1-{{z}^{2}} \right)}^{n}}$$
The fact that this is indeed true can be shown quite easily by showing that this representation satisfies the DE (as we have done), the orthogonality condition (which we haven’t touched yet) and is indeed a polynomial (which is obvious). From the theory of analytic functions we know
$${{f}^{\left( n \right)}}\left( z \right)=\frac{n!}{2\pi i}\int\limits_{C}^{{}}{\frac{f\left( s \right)}{{{\left( s-z \right)}^{n+1}}}ds}$$
Hence
$${{P}_{n}}\left( z \right)=\frac{1}{2\pi i}\int\limits_{C}^{{}}{\frac{{{\left( {{s}^{2}}-1 \right)}^{n}}}{{{2}^{n}}{{\left( s-z \right)}^{n+1}}}ds}$$
where C is a contour encircling the point z once. The validity of our results linking back to the DE requires $\left| z \right|\le 1$ and so for convenience we choose the contour C to be $\left| s \right|=1$ .
Next let’s generate the following function
$$g\left( t,z \right)=\sum\limits_{n=0}^{\infty }{{{P}_{n}}\left( z \right){{t}^{n}}}=\frac{1}{2\pi i}\int\limits_{\left| s \right|=1}^{{}}{\sum\limits_{n=0}^{\infty }{\frac{{{t}^{n}}}{{{2}^{n}}}}\frac{{{\left( {{s}^{2}}-1 \right)}^{n}}}{{{\left( s-z \right)}^{n+1}}}ds}=\frac{1}{2\pi i}\int\limits_{\left| s \right|=1}^{{}}{\frac{2}{t+2s-t{{s}^{2}}-2z}ds}
\\=-\frac{2}{t}\frac{1}{2\pi i}\int\limits_{\left| s \right|=1}^{{}}{\frac{1}{{{s}^{2}}-\frac{2}{t}s+\frac{2z}{t}-1}ds}$$
Note the zeros of the polynomial are
$${{s}^{2}}-\frac{2}{t}s+\frac{2z}{t}-1=\left( s-{{s}_{+}} \right)\left( s-{{s}_{-}} \right)\Rightarrow {{s}_{\pm }}=\frac{1\pm \sqrt{1-2zt+{{t}^{2}}}}{t}$$
and we need to find if any of these lie within the contour. Consider
$$q\left( s \right)=\underbrace{t{{s}^{2}}+2z-t}_{f\left( s \right)}\underbrace{-2s}_{g\left( s \right)}=f\left( s \right)+g\left( s \right)$$
We have along the contour $\left| s \right|=1$
$$\left| g\left( s \right) \right|=\left| -2s \right|\underset{\left| s \right|=1}{\mathop{=}}\,2\underset{\left| z \right|<1}{\mathop{\ge }}\,2\left|z\right|\ge \left| f\left( s \right) \right|\underset{\left| s \right|=1}{\mathop{=}}\,\left| t{{s}^{2}}+2z-t \right|$$
By Rouche’s theorem $q\left( s \right)$ has the same number of zeros inside $\left| s \right|<1$ as does $g\left( s \right)$ which has only one. Therefore only one zero will contribute to the integral and it will be the one with the inferior magnitude, i.e. ${{s}_{-}}$ . Calculating the integral via the residue at ${{s}_{-}}$ therefore yields
$$g\left( t,z \right)=\sum\limits_{n=0}^{\infty }{{{P}_{n}}\left( z \right){{t}^{n}}}=-\frac{2}{t\left( {{s}_{-}}-{{s}_{+}} \right)}=\frac{1}{\sqrt{1-2zt+{{t}^{2}}}}$$
Best Answer
Keep in mind that your generating function is a function of two variables, so when you are taking partial derivatives with respect to $x$ and $t$ and they have different effects.
In fact, you can see that if $f(x,t)=(1-2xt+t^2)^{-1/2}$ then $$ f_x(x,t)=\frac{\partial f(x,t)}{\partial x}=\sum_{n=0}^{\infty}t^nP_n'(x) $$ whereas $$ f_t(x,t)=\frac{\partial f(x,t)}{\partial t}=\sum_{n=0}^{\infty}nt^{n-1}P_n(x). $$ Thus, the recurrence relation is equivalent to the following derivative relationship: $$ \frac{f_x(x,t)-f_x(x,0)}{t}-tf_x(x,t)=2tf_t(x,t)+f(x,t)\qquad (\star) $$ since if you substitute in the power series, the $t^n$ coefficient on the left side equals $P'_{n+1}(x)-P'_{n-1}(x)$ and on the right side, the $t^n$ coefficient equals $(2n+1)P_n(x)$.
The equation $(\star)$ is a partial differential equation for $f(x,t)$, so all you need to do is compute the two derivatives $f_x(x,t)$ and $f_t(x,t)$, and verify that both sides of $(\star)$ match.