Creo que su inicio es correcta. Considere la posibilidad de
$$
F(a_0, \dots, a_N) = \int \limits_{-1}^1 \left( \sum \limits_{k=0}^N a_k x^{2k+1} - \sin(\pi x) \right)^2 dx
$$
El uso de la $L^2$ producto escalar
$$
\langle f, g \rangle := \int \limits_{-1}^1 f(x) g(x) dx
$$
y el correspondiente $L^2$norma $\| f \|_2 := \sqrt{\langle f, f \rangle}$, puede volver a escribir a partir de la expresión como (utilizando la abreviatura $\sum \dots$ para la suma y suponiendo que todos los coeficientes son reales, es decir, $a_j \in \mathbb{R}$ todos los $j$)
$$
F(a_0, \dots, a_N) = \int \limits_{-1}^1 \left( \sum \limits_{k=0}^N a_k x^{2k+1} - \sin(\pi x) \right)^2 dx = \left\| \sum \dots - \sin(\pi x) \right\|_2^2 = \\
= \langle \sum \dots \sum \dots \rangle - 2 \langle \sum \dots, \sin(\pi x) \rangle + \langle \sin(\pi x), \sin(\pi x) \rangle
$$
Debido a la linealidad del producto escalar, podemos escribir esto como
$$
\sum \limits_{k,j=0}^N a_k a_j \int \limits_{-1}^1 x^{2(k+j)+2} dx - 2 \sum \limits_{j=0}^N a_j \int \limits_{-1}^1 \sin(\pi x) x^{2j+1} dx + \int \limits_{-1}^1 \sin^2(\pi x) dx
$$
en el que se evalúa a
$$
\sum \limits_{k,j=0}^N a_k a_j \int \limits_{-1}^1 x^{2(k+j)+2} dx - 2 \sum \limits_{j=0}^N a_j \int \limits_{-1}^1 \sin(\pi x) x^{2j+1} dx + 1
$$
La primera integral, simplemente tiene el valor de $\frac{2}{2(k+j)+3}$ por la integración de primaria. Ahora factorización de la suma con el índice de $j$ finalmente se obtiene:
$$
2 \sum \limits_{j = 0}^N a_j \left( \sum \limits_{k=0}^N \frac{a_k}{2(k+j)+3} - \int \limits_{-1}^1 \sin(\pi x) x^{2j+1} dx \right) + 1
$$
La expresión entre corchetes de los rendimientos exactamente el término $I_j$ que estás buscando.