**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](https://i.stack.imgur.com/IfLEm.jpg) and $f\_3$ to $f\_{15}$: [f\_3 to f\_15](https://i.stack.imgur.com/de5pz.jpg), 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$.



Orthogonal system of functions ordered by norm of second derivative
Here is the proof$^\ast$ (and answer to my question) which I promised.


The functions $f\_1$ and $f\_2$ form the set of affine functions and we can define
\begin{align}
 f\_1(x) = \frac{1}{\sqrt{2}}, \quad f\_2(x) = \sqrt{\frac{3}{2}} x.
\end{align}
All subsequent functions are given (in correct order and normalized) by
\begin{align}
\sqrt{1+\frac{\cos(\nu\_{1+2i})^2}{\cosh(\nu\_{1+2i})^2}} f\_{1+2i} & = \cos(\nu\_{1+2i}x)+\frac{\cos(\nu\_{1+2i})}{\cosh(\nu\_{1+2i})} \cosh(\nu\_{1+2i}x) \\
 \sqrt{1-\frac{\sin(\nu\_{2+2i})^2}{\sinh(\nu\_{2+2i})^2}} f\_{2+2i} & = \sin(\nu\_{2+2i}x)+\frac{\sin(\nu\_{2+2i})}{\sinh(\nu\_{2+2i})} \sinh(\nu\_{2+2i}x), \quad i \in \mathbb{N}
\end{align}
where $\nu\_{1+2i}$ is the $i$-th (positive) root of $\cos(\nu)\sinh(\nu)+\cosh(\nu)\sin(\nu)$ and $\nu\_{2+2i}$ is the $i$-th (positive) root of $\cos(\nu)\sinh(\nu)-\cosh(\nu)\sin(\nu)$. Further more, we have
\begin{align}
 \int\_{-1}^1 f\_i'' f\_j'' = \delta\_{ij} \nu\_i^4.
\end{align}


$^\ast$: At one point we assume $f\_k$ to be sufficiently smooth. I guess this follows easily within a standard argument using calculus of variations, but every derivation I found so far, just makes the same assumption.


**Proof:**


Using Lagrange Multipliers (this part is thanks to Dirk) in order to obtain $f\_k$, $k > 3$, we form the Lagrangian
\begin{align}
 \mathcal{L}^{(k)}(g,\lambda^{(k)},\mu^{(k)}) & = \int\_{-1}^1 |g''|^2 + \sum\_{i = 1}^{k-1} \lambda\_i^{(k)} \int\_{-1}^1 g f\_i 
 + \mu^{(k)} \left( \int\_{-1}^1 |g|^2 - 1 \right).
\end{align}
The derivatives with respect to $\lambda\_1,\ldots,\lambda\_{k-1}$ give the orthogonality conditions, 
whereas $\mu$ gives the normalization. Variation of the function $g$ provides 
\begin{align}
 2 \int\_{-1}^1 g'' h'' + \sum\_{i = 1}^{k-1} \lambda\_i^{(k)} \int\_{-1}^1 f\_i h 
 + 2 \mu^{(k)} \int\_{-1}^1 g h = 0.
\end{align}
Using partial integration twice (at this points we actually assume $f$ to be sufficiently smooth)
we get
\begin{align}
 & 2 \int\_{-1}^1 g^{(4)} h + 2 g^{(2)} h|\_{-1}^1 - 2 g^{(3)} h|\_{-1}^1 \\
 & + \sum\_{i = 1}^{k-1} \lambda\_i^{(k)} \int\_{-1}^1 f\_i h 
 + 2 \mu^{(k)} \int\_{-1}^1 g h = 0.
\end{align}
Since $h \in C^\infty$ is arbitrary, we get the pointwise identity
\begin{align}
 \mathcal{P}^{(k)}(g) := 2g^{(4)} + \sum\_{i = 1}^{k-1} \lambda\_i^{(k)}f\_i
 + 2 \mu^{(k)} g = 0
\end{align}
as well as the boundary conditions
\begin{align}
 g^{(2)}(\pm 1) = 0, \quad g^{(3)}(\pm 1) = 0.
\end{align}
We inductively proof now that $\lambda^{(k)} \equiv 0$, $k \geq 3$. For $j < k$, by orthonormality conditions,
we can simplify
\begin{align}
 0 = \int\_{-1}^1 \mathcal{P}^{(k)}(g) f\_j & = 
 2 \int\_{-1}^1 g^{(4)} f\_j + \sum\_{i = 1}^{k-1} \lambda\_i^{(k)} \int\_{-1}^1 f\_i f\_j 
 + 2 \int\_{-1}^1 \mu^{(k)} g f\_j \\
 & = 
 2 \int\_{-1}^1 g^{(4)} f\_j + \lambda\_j^{(k)}
\end{align}
Using the boundary conditions for both $f$ and $g$, this is equivalent to
\begin{align}
 0 = \int\_{-1}^1 g f\_j^{(4)} + \lambda\_j^{(k)}
\end{align}
In the case of $k = 3$, both $f\_1^{(4)} \equiv f\_2^{(4)} \equiv 0$. 
Hence, $\lambda\_j^{(3)} = 0$, $j = 1,2$. We remain with
\begin{align}
 g^{(4)} + \mu^{(k)} g = 0
\end{align}
In the $k$-th inductive step, we can hence assume that $f\_j^{(4)} = -\mu^{(j)} f\_j$.
Again, by orthogonality, we conclude $\lambda\_j^{(k)} = 0$ for (all) $k$. Thereby,
$f\_k$ is the one normalized solution of the simple differential equation 
\begin{align}
 f\_k^{(4)} + \mu^{(k)} f\_k = 0, \quad f\_k^{(2)}(\pm 1) = 0, \ f\_k^{(3)}(\pm 1) = 0.
\end{align}
Searching for its solutions, $f\_k$ can be assumed to be either symmetric or skew-symmetric. 
If neither was the case, then due to the symmetry of the problem, $g(x) := f\_k(-x)$ fulfills the same differential equation
(for the same $\mu^{(k)}$). Hence, $g + f\_k$ and $g - f\_k$ are 
symmetric and skew-symmetric, respectively, fulfill the differential equation and span the same space
as $g$ and $f\_k$. For simplicity, we substitute $\nu = \sqrt[4]{-\mu^{(k)}}$. 
The unscaled symmetric and skew-symmetric solutions are given by
\begin{align}
 \cos(\nu x) + a \cosh(\nu x), \quad \sin(\nu x) + b \sinh(\nu x),
\end{align}
respectively.
In order to satisfy the boundary conditions,
one calculates that
\begin{align}
 \begin{pmatrix}
 -\cos(\nu) & \cosh(\nu) \\ \sin(\nu) & \sinh(\nu)
 \end{pmatrix}
 \begin{pmatrix}
 1 \\ a
 \end{pmatrix}
 = 0, \quad
 \begin{pmatrix}
 -\sin(\nu) & \sinh(\nu) \\ -\cos(\nu) & \cosh(\nu)
 \end{pmatrix}
 \begin{pmatrix}
 1 \\ b
 \end{pmatrix}
 = 0.
\end{align}
The points at which either matrix becomes singular, are the roots $r\_c$ and $r\_s$ of
\begin{align}
t\_c := \cos(\nu)\sinh(\nu)+\cosh(\nu)\sin(\nu), \quad t\_s := \cos(\nu)\sinh(\nu)-\cosh(\nu)\sin(\nu),
\end{align}
respectively. Both have infinitely many (single) roots which can be approximated by $-\pi/4+i\pi$ and $\pi/4+i\pi$, respectively, since
for $k > 7$, these equal the true $i$-th roots $r\_c^{(i)}$ and $r\_s^{(i)}$ up to $20$ decimal places. 
Earlier values are easily found numerically. 
We therefor define $\nu\_{1,2} = 0, \ \nu\_{1+2i} = r^{(i)}\_c$ and $\nu\_{2+2i} = r^{(i)}\_s$.
After some calculus, we obtain
\begin{align}
\sqrt{1+\frac{\cos(\nu\_{1+2i})^2}{\cosh(\nu\_{1+2i})^2}} f\_{1+2i} & = \cos(\nu\_{1+2i}x)+\frac{\cos(\nu\_{1+2i})}{\cosh(\nu\_{1+2i})} \cosh(\nu\_{1+2i}x) \\
 \sqrt{1-\frac{\sin(\nu\_{2+2i})^2}{\sinh(\nu\_{2+2i})^2}} f\_{2+2i} & = \sin(\nu\_{2+2i}x)+\frac{\sin(\nu\_{2+2i})}{\sinh(\nu\_{2+2i})} \sinh(\nu\_{2+2i}x), \quad i \in \mathbb{N}
\end{align}


All these roots are distinct. The orthogonality conditions are easy to verify, considering
 \begin{align}
 \int\_{-1}^1 f\_i f\_j = \frac{1}{{\nu\_i}^4} \int\_{-1}^1 f\_i^{(4)} f\_j = -\frac{1}{{\nu\_i}^4} \int\_{-1}^1 f\_i f\_j^{(4)}
 = \frac{\nu\_j^4}{{\nu\_i}^4} \int\_{-1}^1 f\_i f\_j
\end{align}
since either $\nu\_j \neq \nu\_i$ or $\nu\_j = 0$, $j = 1,2$, for $j < i$. 
One can further calculate
\begin{align}
 \int\_{-1}^1 f\_i'' f\_j'' = \delta\_{ij} \nu\_i^4.
\end{align}
The $f\_k$ are hence in correct order and their second derivates are orthogonal as well.


