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.
In this answer, the coefficients are computed analytically.
We have $\sin\pi x=\sum_{n=1}^\infty (4n-1)a_n P_{2n-1}(x)$ with $a_n=\int_0^1 P_{2n-1}(x)\sin\pi x\,dx$.
Using $(2n+1)P_n(x)=\big[P_{n+1}(x)-P_{n-1}(x)\big]'$ and integration by parts, we get
$$\left.\begin{aligned} C_n&:=\int_{-1}^1 P_n(x)\cos\pi x\,dx \\S_n&:=\int_{-1}^1 P_n(x)\sin\pi x\,dx \end{aligned}\right\} \implies \left\{\begin{aligned} (2n+1)C_n&=\pi(S_{n+1}-S_{n-1}) \\(2n+1)S_n&=\pi(C_{n-1}-C_{n+1}) \end{aligned}\right.$$
which gives a computation of $a_n=S_{2n-1}/2$, with $C_0=C_1=S_0=0$ and $S_1=2/\pi$.
Alternatively, we have the following expression: $$a_n=\sum_{k=0}^{n-1}\frac{(-1)^k}{(2k)!}\frac{(2n+2k-1)!}{(2n-2k-1)!}\frac{1}{2^{2k}\pi^{2k+1}}.$$ A (long) way to get it: take $a_n=\frac12\int_{-1}^1$ and use Rodrigues' formula for $P_{2n-1}(x)$, then integrate by parts $2n-1$ times, then use Poisson's integral for Bessel functions to get $$a_n=\frac{(-1)^{n-1}}{\sqrt2}J_{2n-\frac12}(\pi)=(-1)^{n-1}j_{2n-1}(\pi),$$ and finally the explicit expression for the spherical Bessel function.