In what ways is the standard Fourier basis optimal? Is there any convincing sense in which the standard trigonometric basis for
the space $V$ of square-integrable real-valued functions on $[-\pi,\pi]$ is
optimal among all the orthonormal bases?
(If this question strikes you as overly subjective, replace "convincing" by "described in the existing mathematical literature".)
Here's an example of the kind of thing that might be true (but please don't
confuse this concrete question with the deliberately vague question that
I'm really interested in!).
Consider the sequence of functions 
$\frac{1}{\sqrt{2}}$, $\cos t$, $\sin t$, $\cos 2t$, $\sin 2t$, $\dots$.
This is an orthonormal basis for $V$ with respect to the inner product
$\langle f,g \rangle = \frac{1}{\pi} \int_{-\pi}^{\pi} f(t)g(t) \: dt$.
The functions have respective total variation 0,4,4,8,8,$\dots$.
I wishfully conjecture (on the basis of no evidence whatsoever)
that for any finite set $\{f_1,\dots,f_n\}$ of orthonormal elements 
of $V$, the total variation of $f_i$, summed over $i$ going from 1 to $n$,
is at least $4 \lfloor n^2/4 \rfloor$
(which is the sum of the first $n$ terms of the sequence $0,4,4,8,8,\dots$).
In this sense, the standard Fourier basis is a "least wiggly" basis.
 A: As you note your question is vague.
There are certainly several ways in which the standard trigonometric system is not optimal among all orthonormal bases (ONB). Let me describe the general philosophy behind these results. Given a set of coefficients ${a_n}$ one would like to understand the function 
$f(x) := \sum a_n \phi_n(x)$
One might start by trying to understand the distribution function of $f$. Generally, the best distributional estimates one can hope for would be those that are implied when $\phi_n$ are independent random variables (the reason this is the best you can hope for is that one could choose the coefficients randomly, from which one can derive lower bounds on the distribution function similar to the case of independent random variables). Of course, most of the time this is far too optimistic since a system of independent random variables is far from complete (for instances, the inequality $||f||_{3} \ll ||f||_{2}$ holds for every function in the span of a system of independent random variables, but certainly doesn't hold in the span of any complete ONB!).  None-the-less, even if one is interested in complete ONB's, there are various natural distributional inequalities of similar strength to the case of independent random variables, after accounting for the obvious ways in which they might fail (and which have a range of interesting consequences), that do not hold for the trigonometric basis, but can be shown to hold for other (and indeed generic) complete ONBs. 
[I have a paper on the topic of constructing orthonormal basis that behave better than the trigonometric system in certain ways: http://arxiv.org/pdf/1205.2420v1.pdf (see Proposition 3 there for an explicit statement of some of the distributional inequalities I refer to above). ]
[I removed a construction of a potential counterexample to your second question, after realizing it wasn't a counterexample]
A: Here's something similar to what you conjecture.  Let's work instead with functions $\mathbb{R} / 2\pi \mathbb{Z} \to \mathbb{C}$; one can do similar things in the real-valued case but I find the complex case simpler to write down. 
Say we want to construct an ONB $f_0, f_1, \ldots$ of $L^2(\mathbb{R} / 2\pi \mathbb{Z})$ consisting of smooth functions that minimize $\| f_n' \|_2$ at each step.  (So in particular, this will minimize $\sum_{k=0}^n \| f_k'\|_2^2$ and other such functionals similar to what you proposed looking at.) Clearly we need $f_0$ to be constant.  Then, expand $f_1$ as a Fourier series:
$$
f_1(x) = \sum_{k\in \mathbb{Z}} a_k e^{ikx}.
$$
We must have $a_0 = 0$ for $f_1$ to be orthogonal to the constant function $f_0$.  To normalize $f_1$ we need 
$$
\sum_{k \in \mathbb{Z},\ k \neq 0} |a_k|^2 = 1
$$ 
(with an appropriate normalization of the norm), and since $f_1$ is smooth,
$$
\| f_1' \|_2^2 = \sum_{k\neq 0} k^2 |a_k|^2.
$$
Clearly, to minimize $\| f_1'\|_2$ subject to $\|f_1\|_2 = 1$, we need $a_k = 0$ for $|k| \ge 2$, so $f_1$ is a linear combination of $e^{ix}$ and $e^{-ix}$.  We'll next want $f_2$ to be an orthogonal linear combination of $e^{ix}$ and $e^{-ix}$.  At the next step, $f_3$ will only have nonzero Fourier coefficients for $|k| \le 2$, and minimizing $\| f_3 \|_2$ forces $f_3$ to be a linear combination of $e^{2ix}$ and $e^{-2ix}$.
Continuing in this way we get that each $f_n$ will be a linear combination of $e^{\lceil n/2 \rceil i x}$ and $e^{- \lceil n/2 \rceil i x}$. We see in particular that $\| f_n' \|_2^2 = \lceil n/2 \rceil$ for each $n$ is the best one can do, and is almost uniquely achieved by the exponential basis.
A: Maybe the question is a bit vague so if this doesn't address whatever was intended, it could be forgivable.
The heat equation is
$$
\frac{\partial u}{\partial t} = c \frac{\partial^2 u}{\partial x^2}
$$
(and then you can add partial derivatives with respect to additional space dimensions if there are any).  Suppose $u(x,t)$ is defined for $0\le x\le 2\pi$ and $t\ge 0$, and the initial temperature is $f(x)$.  Write $f(x)$ as an infinite sum of sine and cosine functions of $x$.  Then each sine or cosine term multiplied by a suitable exponential function of $t$ solves the heat equation, and the rate of decay of the exponential function depends on the frequency of the sine or cosine function.  Since each of these is a solution of the heat equation, so is the infinite linear combination of them.
If you want $u(x,t)=a(x)b(t)$ to be a solution of the heat equation, then $a(x)$ has to be a trigonometric function of $x$ and $b(t)$ has to be an exponential function of $t$.
I think that's what Fourier had in mind.
