# Problem with Legendre-Fourier series for sinx when the number of terms approaches infinity

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.

• Obviously, this is a numerical error. It's not a good idea to involve factorials of numbers greater than 20. Dec 29, 2021 at 7:36
• @SergeiLytkin Oh... now I realized how big 20! is... so this is due to unstability of the computation of large numbers arising from numerical errors. Do I correctly understand what you mean?
– t0d4
Dec 29, 2021 at 8:48
• Yep. There is a more reasonable way to compute the values of the Legendre polynomials via the recurrence relation $(n+1)P_{n+1}(x) = (2n+1)xP_n(x) - n P_{n-1}(x)$. Dec 29, 2021 at 9:09
• @SergeiLytkin Thank you so much for pointing this out. I updated the notebook so that it uses the recurrence relation to calculate a Legendre polynomial of order n. However, the plot still seems to be wrong... Is there a factor that makes this calculation inherently impossble(such as an overflow)? Sorry to bother you so much.
– t0d4
Dec 29, 2021 at 12:36

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.

• Thank you so much for showing me that this problem is analytically solveable. Although I'm not familiar enough with mathematics to fully understand this answer (my major is not mathematics...), this is really interesting.
– t0d4
Feb 1 at 12:54