Orthogonal system of functions ordered by norm of second derivative Problem setting: 
Let $\Omega = [-1,1] \subset \mathbb{R}$ be an interval and consider the space of infinitely differentiable functions, that is $C^{\infty}$. 
We successively define the sequence $f_k \in C^{\infty}, \ k \in \mathbb{N}$, by
\begin{align}
\|f_k''\|_{L_2(\Omega)} & = \inf_{g \ \in \ C^{\infty}}\|g''\|_{L_2(\Omega)} \\
& \mbox{subject to:} \ g \perp \ \mathrm{span}(f_1,\ldots,f_{k-1}), \ \|g\|_{L_2} = 1
\end{align}
for the standard $L_2$ scalar product, such that
\begin{align}
\int_{\Omega} f_i \ f_j = \delta_{ij}.
\end{align}
Question: Is there a numerically feasible description of these functions, as well as their first and second derivatives? And what is the space spanned by these functions?
Comments on the problem: Quite obviously, $f_1$ and $f_2$ form the set of affine functions and I guess all further functions are unique. 
By numerically feasible description I mean sorts of functions which are (numerically) evaluable with reasonable computational complexity at any one point (e.g. polynomials, trigonometrical functions, recursively defined functions, ...). Any further theoretical insight on these functions is helpful nonetheless. 
If one replaces $\|\cdot''\|$ by $\|\cdot'\|$, then the $f_k$ are seemingly (ignoring their order and modulo scaling) given by the system
\begin{align}
1, \quad\sin((2i + 1)\pi/2 \ x), \quad \cos((2i + 2)\pi/2 \ x), \quad i \in \mathbb{N}.
\end{align}
All these functions have derivative $0$ at the interval border. 
For the original problem, by an SVD of a discretized version of the problem, one can get a notion of how these functions look like (the plots are pretty exact). Here are functions $f_3$ to $f_7$: f_3 to f_7 and $f_3$ to $f_{15}$: f_3 to f_15, from which one may observe several properties. For example, the second derivative at the border seems to be $0$.
Regarding the space $C^{\infty}$, I think it is unlikely, yet I am not sure, that the original phrasing of the problem is well-defined. In case of need, one may replace it by $C^2$ or even $H^2$ (the space of twice weakly differentiable functions). Note that the latter one is embedded into $C^1$.
 A: The optimization problem for $f_3$ is
$$
\min_g \int_{-1}^1 |g''(x)|^2dx \quad \text{s.t.}\quad \int_{-1}^1 g(x)dx = 0,\ \int_{-1}^1 x g(x)dx = 0,\ \int_{-1}^1 |g(x)|^2dx = 1.
$$
Using Lagrange Multipliers to solve this optimization problem you form the Lagrangian
$$
\mathcal{L}(g,\lambda_1,\lambda_2,\mu) = \int_{-1}^1 |g''(x)|^2dx + \lambda_1 \int_{-1}^1 g(x)dx + \lambda_2\int_{-1}^1 xg(x)dx + \mu\left(\int_{-1}^1|g(x)|^2dx-1\right).
$$
Now derive $\mathcal{L}$ with respect to all variables and set them to zero to get the necessary (and, due to convexity of the problem, also sufficient) conditions for optimality.
The derivatives for $\lambda_1,\lambda_2$ and $\mu$ just give the constraints and for the derivative for the "function variable" just do a formal calculation to calculate the directional derivative
$$
\lim_{t\to 0} \frac{\mathcal{L}(g+th,\dots) - \mathcal{L}(g,\dots)}{t}.
$$
Using integration by parts this shows that $f_3$ has to solve the boundary value problem
$$
g^{(4)}(x) = -\lambda_1 - \lambda_2 x - 2\mu g(x),\quad g^{(2)}(\pm 1) = 0,\quad g^{(3)}(\pm 1) = 0 
$$
for some Lagrange multiplies $\lambda_{1/2},\mu$ and also fulfill the additional constraints. There is a closed form solution for this problem, but I did not push the calculation through…
Proceeding similarly gives $f_k$ as solution of
$$
g^{(4)}(x) = -\lambda_1 - \lambda_2 x - \lambda_3 f_3(x) - \cdots - \lambda_{k-1}f_{k-1}(x) - 2\mu g(x),\quad g^{(2)}(\pm 1) = 0,\quad g^{(3)}(\pm 1) = 0 
$$
for some Lagrange multiplies $\lambda_{l},\mu$, $l=1,\dots,k-1$. While it seem not too hard to write down $f_3$ explicitly, $f_4$ seems to be quite a mess.
