Note that $a_0$ doesn't really have a separate formula:
If you plug in $n=0$ into $a_n$ you get the same integral.
For an intuition as to why these formula are what they are, say you wanted to calculate a Fourier series for $\cos(nx$). You know that all the coefficients except for $b_n$ must be zero. So how do you make that happen? You use the fact that:
$$\int_{-\pi}^\pi \cos(nx)\cos(mx)dx=0 \ \ \bigg| \ \ m\neq n$$
$$\frac{1}{\pi}\int_{-\pi}^\pi \cos(nx)\cos(mx)dx=1 \ \ \bigg| \ \ m= n=0$$
Doing the same thing for $\sin$, you'll see that only way you can get the coefficients to be consistent with $f(x)=\cos(nx)$ or $f(x)=\sin(nx)$ is to choose the formula's you brought up.
This answer is mostly for students who used an algebra approach. I don't know if Fourier himself thought up the series this way, but it is common today. I left a lot of steps out and mainly showed ideas that I struggled with when I first tried to motivate the Fourier Series.
I'll start off by observing a trigonometric polynomial:
$T(x) = c_0 + c_1 \cos(x) + c_2 \cos(2x)+...+c_n \cos(n x) + d_1 \sin(x) + ... + d_n \sin(n x)$
where $c_n$ and $d_n$ are some non-zero value.
The goal is write the orthogonal basis, from there I can find the coefficients.
So, I need declare the inner product:
$<\textbf{f},\textbf{g}> = \int_0^{2\pi} f(x)g(x)dx$
Before I can obtain a orthogonal basis, I should first get the orthonormal basis by using the Gram–Schmidt process. Where $\|\textbf{g}_0\| = \|\textbf{g}_1\| =\|\textbf{g}_2\| = ... =\|\textbf{g}_n\| = 1$.
$\|\textbf{f}\|^2 = <\textbf{f}><\textbf{f}> = 2\pi$
Thus, $\|\textbf{f}\| = \sqrt{2\pi}$
and
$e_0 = \frac{\textbf{f}}{\|\textbf{f}\|}$
$\textbf{g}_0 = e_0 = \frac{1}{\sqrt{2\pi}}$
$\textbf{g}_1 = e_1 = \frac{1}{\sqrt{\pi}}\cos(x)$
...
$\textbf{g}_n = e_n = \frac{1}{\sqrt{\pi}}\cos(nx)$...
$\textbf{g}_{n+1} = e_{n+1} = \frac{1}{\sqrt{\pi}}\sin((n+1)x)$...
$\textbf{g}_{2n} = e_{2n} = \frac{1}{\sqrt{\pi}}\sin(nx)$
The orthogonal basis yields:
$a_0 = \frac{2}{\sqrt{2\pi}}<\textbf{f},\textbf{g}_0>$ (Use 2 because it makes generalizing the coefficients possible.)
$a_1 = \frac{1}{\sqrt{\pi}}<\textbf{f},\textbf{g}_1>$...
$a_n = \frac{1}{\sqrt{\pi}}<\textbf{f},\textbf{g}_n>$...
$b_n = \frac{1}{\sqrt{\pi}}<\textbf{f},\textbf{g}_{2n}>$
Now that everything is divided into sines and cosines I can get the coefficients:
$a_n = \frac{1}{\pi}\int_0^{2\pi} f(x)\cos(nx)dx $
and
$b_n = \frac{1}{\pi}\int_0^{2\pi} f(x)\sin(nx)dx $
Best Answer
Starting with $$f(x) = \frac{a_0}{2}+\sum_{n=1}^\infty [a_n \cos(nx) + b_n \sin(nx)],$$ to find the $b_n$, multiply both sides of the equation by $\sin(mx)$ for an arbitrary but fixed positive integer $m$, and integrate from $x=0$ to $x=\pi$. Using the fact that $$\int_0^\pi \sin(nx)\sin(mx)\,dx=\begin{cases}0, &n\not=m,\\ {\pi\over 2}, &n=m,\end{cases}$$ you get $b_m={2\over \pi}\int_0^\pi f(x)\sin(mx)\,dx$ for each $m=1,2,\dots$ Since $m$ was arbitrary, you can change this to an $n$ and get the formula for the coefficients for the sine terms.
A similar argument for the cosine terms establishes the formula for the $a_n$. Here, use the fact that $$\int_0^\pi \cos(nx)\cos(mx)\,dx=\begin{cases} 0, &n\not=m,\\ {\pi\over 2}, &n=m,\end{cases}$$ which leads to $$a_n={2\over \pi}\int_0^\pi f(x)\cos(nx)\,dx, \quad n=0,1,2,\dots$$
Note: The only purpose of the factor ${1\over 2}$ on the $a_0$ term is so that the formula for $a_0$ will match the pattern of the formula for $n=1,2,\dots$