# How can I make FourierCoefficient et al. handle the zero-order term correctly?

Say that I would like to use Mathematica to step through one of the most famous solutions to the Basel Problem, using Parseval's Identity,

$$\sum_{n = -\infty}^\infty \left|c_n\right|^2 = \frac{1}{2\pi} \int_{-\pi}^{\pi}dx\,\left|\,f(x)\right|^2$$

where $c_n$ is the $n$th (complex) Fourier coefficient of $f(x)$. I compute this with FourierCoefficient for $f(x) = x$:

In[1]:= FourierCoefficient[x, x, n]
Out[1]= (I (-1)^n)/n


That's correct for any $n \in \mathbb{Z}$ except $0$. If you ask for the specific coefficient, it comes back with the right answer:

In[2]:= FourierCoefficient[x, x, 0]
Out[2]= 0


However, I see no way of getting FourierCoefficient to return a result for general $n$ that correctly handles $0$. The obvious-seeming option of GenerateConditions doesn't do the trick at all:

In[3]:= FourierCoefficient[x, x, n, GenerateConditions -> True]
Out[3]= (I*(-1)^n)/n


This doesn't do you much good if you're trying to carry out that sum, though:

In[4]:= Sum[Abs[FourierCoefficient[x, x, n]]^2, {n, -Infinity, Infinity}]

Power::infy: Infinite expression 1/0^2 encountered. >>

Out[4]= Sum[1/(E^(2*Pi*Im[n])*Abs[n]^2), {n, -Infinity, Infinity}]


The way it's handling the imaginary part of n seems a little eccentric as well, and playing with Assumptions, ComplexExpand and the like did little to resolve that. However, the real problem is the Power::infy message. Now, if you compute the coefficient directly, you get:

In[5]:= term = Integrate[x*Exp[-I*n*x], {x, -Pi, Pi}]/(2*Pi)
Out[5]= (I*(n*Pi*Cos[n*Pi] - Sin[n*Pi]))/(n^2*Pi)


Naïve substitution of n -> 0 still blows up, of course, but it's correct in the limit, but simplification with the assumption that n is an integer runs into the same issue as FourierCoefficient:

In[6]:= Limit[term, n -> 0]
Out[6]= 0

In[7]:= FullSimplify[term, Assumptions -> Element[n, Integers]]
Out[7]= (I*(-1)^n)/n


With the correct set of assumptions and the explicit Limit, you can evaluate the sum:

In[8]:= Assuming[{Element[n, Reals]},
Sum[Limit[Abs[term]^2, n -> k], {k, -Infinity, Infinity}]]
Out[8]= Pi^2/3


This seems like a lot of human effort in order to make Mathematica get the right answer. Is there a way to avoid it?

One way to avoid the singularity in 0-th coefficient is to regularize the original function. The problem is due to it being a polynomial, which is why we get powers of n in denominator by partial integration - the procedure that is obviously failing for the 0-th harmonic. Obviously this is a shortcoming for the current way Mathematica computes Fourier coefficients.

There are many ways one can regularize. For this particular case, I chose the following one:

coef = FourierCoefficient[t*Cos[eps*t], t, n]

(*
-((I (-1)^n n ((eps^2 - n^2) \[Pi] Cos[eps \[Pi]] -
2 eps Sin[eps \[Pi]]))/((eps - n)^2 (eps + n)^2 \[Pi]))
*)


The general procedure then would be to do something like this:

coefsq =
FullSimplify[
#*Conjugate[#] &[ComplexExpand[coef]],
Assumptions ->
Element[n, Integers] && Element[{eps, t}, Reals] && eps > 0
]

(*
(n (-eps + n) (eps + n) \[Pi] Cos[eps \[Pi]] +
2 eps n Sin[eps \[Pi]])^2/((eps^2 - n^2)^4 \[Pi]^2)
*)


And then rely on Mathematica's symbolic capabilities to perform the sum:

sum = Sum[coefsq, {n, -Infinity, Infinity}]

(*
(4 eps^3 \[Pi]^3 + 6 eps \[Pi] Cos[2 eps \[Pi]] - 3 Sin[2 eps \[Pi]]
+ 6 eps^2 \[Pi]^2 Sin[2 eps \[Pi]])/(24 eps^3 \[Pi])
*)


The last step is to take the limit eps -> 0, to remove the regulator:

Limit[sum, eps -> 0]

(* \[Pi]^2/3 *)


which is the correct result.

This procedure isn't perhaps too simple, and also may not work for more complicated functions (at least symbolically), but it at least avoids the main problem you reported - the singularity for n=0.

• +1; I definitely like the regularization idea. – Pillsy Nov 3 '15 at 18:51

Evidently, one answer to this question is to upgrade Mathematica. Version 11.2 returns the following:

FourierCoefficient[x, x, n]
(* Piecewise[{{0, n == 0}}, (I*(-1)^n)/n] *)


It's a nice improvement.