Find the continuous function defined over $[0,1]$ such that $\int_0^1f(x){\rm d}x=1, \int_0^1 xf(x){\rm d}x=a$ and $\int_0^1x^2f(x){\rm d}x=a^2 .$ Find the continuous function defined over $[0,1]$ such that $$\int_0^1 x^kf(x){\rm d}x=a^k,$$ where $k=0,1,2.$
My Thought
Indeed, the conditions imply that $$\int_0^1f(x)=1,~~~\int_0^1 xf(x){\rm d}x=a,~~~\int_0^1 x^2f(x){\rm d}x=a^2.$$ Thus, it's easy to obtain that $$\begin{align*}\int_0^1(1-x)^2f(x){\rm d}x&=\int_0^1(1-2x+x^2)f(x){\rm d}x\\&=\int_0^1f(x){\rm d}x-2\int_0^1xf(x){\rm d}x+\int_0^1 x^2f(x){\rm d}x\\&=1-2a+a^2\\&=(1-a)^2\end{align*}.$$
But how to go on with this? I'm stuck here. Please offer some help. Thanks.
 A: There are infinitely many functions with this property. If you are trying to find one such function take $f(x)$ to be a polynomial of degree 2 and solve for the coefficients from the given conditions.
A: Since
$$
\int_0^1 {f(x)dx}  = 1
$$
then f(x) - if assumed to be non-negative - is a probability distribution function that is defined by its moments.
Your problem goes under the Hausdorff moment problem class.
You do not specify what the moments higher than $2$ shall be, so let's examine two different cases.
a) case $k = 0..\infty$
Let's consider first the case in which the all the moments are specified to be equal to $a^k$.
Then given
$$
\int_0^1 {x^{\,k} f(x)dx}  = a^{\,k} 
$$
it is also true that
$$
{{t^{\,k} } \over {k!}}\int_0^1 {x^{\,k} f(x)dx}  = \int_0^1 {{{t^{\,k} x^{\,k} } \over {k!}}f(x)dx}  = {{t^{\,k} a^{\,k} } \over {k!}}
$$
and that
$$
\int_0^1 {\sum\limits_{0\, \le \,k} {{{t^{\,k} x^{\,k} } \over {k!}}} f(x)dx}  = \int_0^1 {e^{\,t\,x} f(x)dx}  = \sum\limits_{0\, \le \,k} {{{t^{\,k} a^{\,k} } \over {k!}}}  = e^{\,a\,t} 
$$
Since your function is defined in $[0,1]$, then we can write
$$
\int_0^1 {e^{\,t\,x} f(x)dx}  = \int_{ - \infty }^\infty  {e^{\,t\,x} f(x)\left( {H(x) - H(x - 1)} \right)dx}  = e^{\,a\,t} 
$$
where $H(x)$ is the Heaviside function, and the expression is
the Moment Generating function or Laplace Transform
of $f(x)$ restricted to the $[0,1]$ range and null outside.
Without going on in that approach, if we compute instead the central moments
$$
\eqalign{
  & \overline x  = \int_0^1 {xf(x)dx}  = a  \cr 
  & \mu _{\,k}  = \int_0^1 {\left( {x - \overline x } \right)^{\,k} f(x)dx}  = \int_0^1 {\sum\limits_{0\, \le \,j\, \le \,k} {\left( { - 1} \right)^{\,k - j} \left( \matrix{
  k \cr 
  j \cr}  \right)\overline x ^{\,k - j} x^{\,j} } f(x)dx}  =   \cr 
  &  = \sum\limits_{0\, \le \,j\, \le \,k} {\left( { - 1} \right)^{\,k - j} \left( \matrix{
  k \cr 
  j \cr}  \right)\overline x ^{\,k - j} \int_0^1 {x^{\,j} f(x)dx} }  =   \cr 
  &  = \sum\limits_{0\, \le \,j\, \le \,k} {\left( { - 1} \right)^{\,k - j} \left( \matrix{
  k \cr 
  j \cr}  \right)a^{\,k} }  = 0^{\,k} a^{\,k}  \cr} 
$$
we discover that they are all null.   
That means that, provided $0 \le a \le 1$,  the only non-negative $f(x)$ having the required moments
is the Dirac delta function $\delta(x-a)$.
b) case higher moments null
Assuming that the moments higher than $2$ are null, we can always proceed as above.   
But if you are looking to express $f(x)$ as a polynomial, waiving on the requirement of being non-negative, then we can follow a shorter way.
Putting
$$
f(x) = \sum\limits_{0\, \le \,n\, \le \,h} {c_{\,n} x^{\,n} } 
$$
we get
$$
a^{\,k}  = \int_0^1 {x^{\,k} f(x)dx}  = \int_0^1 {\sum\limits_{0\, \le \,n\, \le \,h} {c_{\,n} x^{\,n + k} } dx}  = \sum\limits_{0\, \le \,n\, \le \,h} {{{c_{\,n} } \over {n + k + 1}}} 
$$
that is
$$
\left( {\matrix{
   1 & {1/2} & {1/3} &  \cdots   \cr 
   {1/2} & {1/3} & {1/4} &  \cdots   \cr 
   {1/3} & {1/4} & {1/5} &  \cdots   \cr 
    \vdots  &  \vdots  &  \vdots  &  \ddots   \cr 
 } } \right)\left( {\matrix{
   {c_{\,0} }  \cr 
   {c_{\,1} }  \cr 
   {c_{\,2} }  \cr 
    \vdots   \cr 
 } } \right) = \left( {\matrix{
   1  \cr 
   a  \cr 
   {a^{\,2} }  \cr 
    \vdots   \cr 
 } } \right)
$$
The matrix is invertible, although I could not find an easy expression for its inverse.
For a limited number of factors, we can invert it numerically and for $h=2$ we obtain for instance
$f(x)=(180a^2 - 180a + 30)x^2 + (192a - 180a^2 - 36)x + (30a^2 - 36a + 9)$
A: Let
$$ f(x)=b_2x^2+b_1x+b_0. $$
Then
$$ \int_0^1f(x)dx=\frac13b_2+\frac12b_1+b_0, \int_0^1xf(x)dx=\frac14b_2+\frac13b_1+\frac12b_0$$
and
$$ \int_0^1x^2f(x)dx=\frac15b_2+\frac14b_1+\frac13b_0. $$
Let
$$ \frac13b_2+\frac12b_1+b_0=1,\frac14b_2+\frac14b_1+\frac12b_0=a,\frac15b_2+\frac14b_1+\frac13b_0=a^2 $$
and then
$$  b_0=3 \left(10 a^2-12 a+3\right),b_1=-12 \left(15
   a^2-16 a+3\right),b_2=30 \left(6 a^2-6 a+1\right). $$
So
$$ f(x)=30 \left(6 a^2-6 a+1\right)x^2-12 \left(15
   a^2-16 a+3\right)x+3 \left(10 a^2-12 a+3\right).$$
