Esta suma es una función de la forma $f(x_1,\ldots,x_n)=x^TMx+\lambda^Tx$ $x=[x_1,\ldots,x_n]^T$ y el vector $\lambda\in\mathbb{R}^n$. Si lleva a cabo los cálculos que toma la forma
$$f(x)=\big(\sum_{i=1}^n{x_i}\big)^2-n\sum_{i=1}^n{x_i^2}-\sum_{i=1}^n{(n-2i+1)x_i}.\qquad \qquad (1)$$
El resultado clave es que podemos escribir (1) de manera equivalente como
$$f(x)=-n\sum_{i=1}^n{\Big(x_i-\frac{1}{n}\sum_{j=1}^{n}{x_j}+\frac{n-2i+1}{2n}\Big)^2}+\frac{1}{4n}{\sum_{i=1}^n(n-2i+1)^2} \qquad\qquad (2)$$
que alcanza el máximo cuando
$$x_i=\frac{1}{n}\sum_{j=1}^{n}{x_j}-\frac{n-2i+1}{2n}\qquad \qquad\qquad\qquad\qquad\qquad\qquad(3)$$
Este valor máximo es alcanzable si $0\leq x_i\leq 1$ o, equivalentemente, si
$\frac{n-2i+1}{2n}\leq \frac{1}{n}\sum_{j=1}^{n}{x_j}\leq\frac{n-2i+1}{2n}+ 1$ todos los $i=1,\ldots,n$
es decir, si
$$\frac{n-1}{2n}\leq \frac{1}{n}\sum_{j=1}^{n}{x_j}\leq\frac{n+1}{2n}$$
La selección de $\frac{1}{n}\sum_{j=1}^{n}{x_j}=\frac{n+1}{2n}$ obtenemos de (3) $x_i=i/n$. Ahora bien, si hemos de sustituir (3) en (2) el valor máximo de
$$\frac{1}{4n}{\sum_{i=1}^n(n-2i+1)^2}=\frac{n^2-1}{12}$$
se deriva.
Una observación interesante sobre el problema es que hay un número infinito de puntos que alcanzar este valor máximo, incluso en $[0,1]^n$. Esto es una consecuencia de (3) donde se puede ver que para cada valor de la suma de $\sum_{i=1}^n{x_i}$ en el intervalo de $[(n-1)/2,(n+1)/2]$ diferentes valores de la óptima $x_i$ se obtienen. Sólo tienen que diferir por $1/n$.