Let $\mathcal{P}_n$ denote the space of polynomials of degree up to $n$. Then $$(P,Q):=\int_{-1}^1 P(x)Q(x)\,dx$$ defines an inner product on $\mathcal{P}_n$, and $l(P):=P(a)$ is a linear functional on it.
By the Riesz representation theorem, there exists $Q_a\in\mathcal{P}_n$ such that $l(P)=(P,Q_a)$. The theorem is more commonly used for Hilbert spaces in functional analysis, and there one needs the functional to be continuous. But for finite dimensional inner product spaces all linear functionals are continuous and one can construct the requisite $Q_a$ explicitly by using an orthonormal basis $E_i$. Namely $$Q_a(x):=\sum_i l(E_i)E_i(x)=\sum_i E_i(a)E_i(x).$$
Indeed, taking the inner product with $E_j$ on both sides gives us $(E_j,Q_a)=E_j(a)$ for all $j\leq n$, and then it follows for all $P$ of degree up to $n$ by linearity.
Of course, the real work is swept under the rug in orthonormalizing the powers $x^i$ into polynomials $E_i(x)$ orthonormal in the above inner product. Fortunately, in this case $E_i(x)$ are just the normalized Legendre polynomials, which are known explicitly.