After I learned about Fourier series expansion, I understand orthogonality of trigonometric functions was the key when I calculate the coefficients of Fourier series. As I knew that Legendre Polynomials are also orthogonal with each other, I came up with a following problem:
Calculate Legendre-Fourier series for $f(x) = \sin{x}$
My solution
First, in this context, inner product of arbitrary functions $g, h$ is defined as $\langle g, h \rangle = \int_{-\pi}^{\pi} g(x)h(x) \mathrm{d}x$.
Next, define $Q_{n}(x)$ as an "adjusted" Legendre Polynomial of order $n$, which means $Q_{n}(x)$ satisfies $\langle Q_{m}(x), Q_{n}(x) \rangle = \delta_{mn}$. So $Q_{n}(x)$ can be written as
$$ Q_{n}(x) = \sqrt{\frac{2n+1}{2\pi}} P_{n}\left( \frac{x}{\pi} \right) $$ where $P_{n}(x)$ is the Legendre Polynominal of order $n$.
Then $f(x)$ is expected to be expanded like $$ f(x) = \sum_{n=0}^{\infty} c_n Q_{n}(x) \tag{1}\label{expected_expansion} $$ where $c_{n}$ are the coefficients.
To calculate the coefficient $c_n$, I can utilize the orthogonality of $Q_{n}(x)$ by multiplying both sides of \eqref{expected_expansion} by $Q_{n}(x)$ and integrating them from $-\pi$ to $\pi$, and this derives
$$ c_n = \int_{-\pi}^{\pi} f(x) Q_{n}(x) \mathrm{d}x $$
Thus, Legendre-Fourier series for $f(x) = \sin{x}$ can be expressed as $$ \sum_{n=0}^{\infty} \left( \int_{-\pi}^{\pi} Q_{n}(x) \sin{x} \mathrm{d}x \right) Q_{n}(x) $$ $$\tag*{$\blacksquare$}$$
Issue
In order to confirm my answer, I defined $\widetilde{f}_{N}(x)$ as $$ \widetilde{f}_{N}(x) = \sum_{n=0}^{N} \left( \int_{-\pi}^{\pi} Q_{n}(x) \sin{x} \mathrm{d}x \right) Q_{n} $$ and plotted $\widetilde{f}_{N}(x)$ in Mathematica to see how well $\widetilde{f}_{N}(x)$ approximates $\sin{x}$ over $-\pi \le x \le \pi$ as $N$ increases. However, although $\widetilde{f}_{N}(x)$ seems to fit $\sin{x}$ nicely when $N \le 20$, it suddenly fluctuates and diverges when $N \ge 21$. (Please look at this notebook) and I have no idea why this is happening.
Since $\sin{x}$ cannot be expressed as a finite sum of polynomials, I think $N$ should approach $\infty$ in order for $\widetilde{f}_{N}(x)$ to be equal to $f(x)$. Are there any errors in my solution or is this a kind of a bug of Mathematica?
Any help would be much appreciated. (Sorry for poor English)
Edit:
Thanks to the @Sergei Lytkin's comment, I now understand the computation of factorial of large numbers was causing this problem. However, even after I modified the notebook so that it uses recurrence formula to calculate Legendre Polynomial, the problem wouldn't vanish. Does anyone have other idea? Thanks in advance.