What is $\langle H \rangle$ in the infinite square well in the state $|\psi \rangle = | x_0 \rangle$, where $0 \le x_0 \le L$? In the 1-dimensional infinite square well (with width $L$) the eigenvalues and the corresponding eigenfunctions are
$$
E_n = \frac{\hbar^2 \pi^2 n^2}{2L^2 m},\qquad \psi_n (x) = \sqrt{\frac{2}{L}}\sin\left( \frac{n\pi}{L}x\right), \qquad x \in [0,L]
$$
The general solution of the time-independent Schrödinger equation in the basis of the Hamitonian is
$$
\psi (x) = \sum_{n=1}^\infty c_n\,\sqrt{\frac{2}{L}}\sin\left( \frac{n\pi}{L}x\right)
$$
Question:
If the particle is in the state $\psi_(x) = \delta(x-x_0)$, where $0 \le x_0 \le L$, i.e. the particle is at the position $x =x_0$ at $t = 0$ with certainity, then what is the expectation value of the energy $\langle H \rangle$ in this state ? or even what is the probability of getting each of the eigenvalues $E_n$ upon measurement ?
I don't know if this is the right way but my attempt is as follows:
I. The coefficients $c_n$ of the expansion using Fourier trick are
$$
c_n = \int_0^L \sqrt{\frac{2}{L}}\sin\left( \frac{n\pi}{L}x\right)\delta(x-x_0)\,dx = \sqrt{\frac{2}{L}}\sin\left( \frac{n\pi}{L}x_0\right)
$$
II. It folows that (I'm not sure if this makes sense)
$$
\langle x_0|x_0\rangle = \sum_{n=1}^\infty |c_n|^2 = \delta(x_0 - x_0) = \infty
$$
III. The expectation value of the energy is
$$
\langle H \rangle = \sum_{n=1}^\infty E_n\,|c_n|^2 = \frac{2}{L}\sum_{n=1}^\infty E_n\,\sin^2\left( \frac{n\pi}{L}x_0\right)
$$
or differently
$$
\langle H \rangle = \int_{-\infty}^\infty \delta(x-x_0)(-\frac{\hbar^2}{2m}\frac{d^2}{dx^2})\delta(x-x_0)\, dx = \cdots
$$
I don't know how to check whether these two expressions of $\langle H \rangle$ are equivalent, how does the first converge and how to compute the second ?

ADDED:
The first expression should be more completely written as
$$
\langle H \rangle = \frac{\langle \psi|H|\psi\rangle}{\langle\psi|\psi\rangle} 
 = \frac{\sum_{n=1}^\infty E_n\,|c_n|^2}{\sum_{n=1}^\infty |c_n|^2} = \frac{\infty}{\infty} = \ ?
$$
and I would suggest to get the limit as the limit of the sequence
$$
\langle H \rangle_m = \left( \frac{\sum_{n=1}^m E_n\,|c_n|^2}{\sum_{n=1}^m |c_n|^2}\right)_{m\in\mathbb N}
$$
Result: it will not converge

 A: 
If the particle is in the state $\psi(x)=\delta(x-x_0) \ldots$

Unfortunately, that's a non-starter.  Delta functions are not square-integrable, and therefore do not represent physical states.  The expansion in terms of energy eigenfunctions is only guaranteed to converge if the function being expanded is square integrable, it's not possible to define expected values in non-normalizable states, etc. As a result, none of the expressions you've written are well-defined.

We can explore your question by considering the following family of functions:
$$\psi_\epsilon(x):= \begin{cases}1/\sqrt \epsilon & x_0-\frac{\epsilon}{2}<x<x_0 + \frac{\epsilon}{2} \\ 0 & \text{else}\end{cases}$$
For all positive $\epsilon$ these functions are square-integrable with norm $\Vert \psi_\epsilon(x) \Vert = 1$.  Therefore, we can plug them into your eigenfunction expansion to find
$$c_n = \sqrt{\frac{2}{\epsilon L}}\int_{x_0-\epsilon/2}^{x_0+\epsilon/2} \sin\left(\frac{n\pi x}{L}\right) \ \mathrm dx = -\sqrt{\frac{2L}{n^2\pi^2\epsilon}}\left[\cos\left(\frac{n\pi x_0}{L} + \frac{n\pi \epsilon}{2L}\right) - \cos\left(\frac{n\pi x_0}{L} - \frac{n\pi \epsilon}{2L}\right)\right]$$
$$= \sqrt{\frac{8L}{n^2\pi^2\epsilon}}\sin\left(\frac{n\pi x_0}{L}\right) \sin\left(\frac{n\pi \epsilon}{2L}\right)$$
and so
$$\left<H\right>_{\psi_\epsilon} = \sum_{n=1}^\infty \frac{n^2\pi^2 \hbar^2}{2mL^2} \cdot \frac{8L}{n^2\pi^2 \epsilon}\sin^2\left(\frac{n\pi x_0}{L}\right) \sin^2\left(\frac{n\pi \epsilon}{2L}\right)$$
$$=\frac{4\hbar^2}{\epsilon mL} \sum_{n=1}^\infty \sin^2\left(\frac{n\pi x_0}{L}\right) \sin^2\left(\frac{n\pi \epsilon}{2L}\right) \rightarrow \infty$$
Now, this sum diverges even if $\epsilon$ is not particularly small, i.e. the particle is not necessarily well-localized.  The reason is because the Hamiltonian is proportional to the second derivative of the wavefunction, so a wavefunction which changes very sharply corresponds to a large expected value of the energy.  Loosely speaking, our wavefunctions change infinitely quickly at the points of discontinuity, meaning that the expected value of the energy is not finite$^\dagger$.
A better example would be to use bump functions - smooth functions with compact support of the form
$$\psi_\epsilon(x) := \begin{cases} A_\epsilon e^{-\epsilon^2/(\epsilon^2-[x-x_0]^2)} & x_0-\epsilon < x < x_0+\epsilon \\ 0 & \text{else}\end{cases}$$
where $A_\epsilon$ is a normalization constant.  These functions are continuous and differentiable at every point, and so they don't suffer from the same problems as our naive trial functions. They are, however, quite tedious to work with in this problem.
What we will find is that for each $\epsilon >0$ the expected value of the energy will be finite - but it will grow without bound as $\epsilon\rightarrow 0^+$.  The reason for this is very similar to the naive case - if the bump function is normalized, then shrinking its width implies making it taller so as to enclose the same area.  But this means that the rate at which the wavefunction changes (and by extension, its second derivative) grows without bound as the wavefunction becomes more and more narrow.  As a result, we should expect that $\lim_{\epsilon\rightarrow 0^+} \left<H\right>_{\psi_\epsilon} \rightarrow \infty$.  I'll trudge through the numerics just to demonstrate this.
To illustrate the point, let $L=\pi,m=\hbar=1$, and $x_0= \pi/2$.  If we define $u\equiv (x-\pi/2)/\epsilon$, we find that
$$A_\epsilon \approx 2.74/\sqrt \epsilon$$
which yields coefficients in the sine expansion
$$c_n = \int_0^\pi\sqrt{\frac{2}{\pi}}\sin(nx) \cdot \psi_\epsilon(x)\mathrm dx$$
$$ \approx \sqrt{\frac{2}{\pi}} \cdot \frac{2.74}{\sqrt{\epsilon}}\epsilon \int_{-1}^1 \sin(n\epsilon u + n \pi/2) e^{-1/(1-u^2)}$$
If $n=2k$ is even, $\sin(n\epsilon u + n\pi/2) = (-1)^k \sin(2k\epsilon u)$ and the coefficient vanishes.  If it is odd, then the coefficients don't have a nice closed form, but can be solved for numerically:

We can then calculate the expected energy as a function of $\epsilon$.  This is the result:

So the main lesson here is that sharp changes in the wavefunction correspond to large expected values of $H$. Discontinuous wavefunctions will have $\langle H\rangle= \infty$, and wavefunctions which are well-localized will have large $\langle H \rangle$ which also tends to infinity as the localization gets sharper and sharper.

$^\dagger$There's nothing unphysical about this statement.  Recall the statistical interpretation of the expected value - it is the average value expected after $N$ identical measurements, as $N\rightarrow \infty$.  A non-finite expected value for the energy simply means that if you make a sufficient number of measurements, then your average result will become arbitrarily large rather than settling down to a finite limit.  In particular, it doesn't mean that any individual measurement will yield an infinite result.
A: The first does not converge. $E_n$ is increasing with $n$, while $\sin^2(n\pi x_0/L)$ is oscillatory and positive-definite. I'll leave the proof of its divergence to you, but you shouldn't expect something to converge when its terms, on average, increase with $n$.
The integral also diverges, as is pretty much always the case when you have the same Dirac delta appearing twice in an integral.
Note that $|x_0\rangle$ is not a normalizable wave function, and as such is not a member of the Hilbert space. You shouldn't expect to get reasonable results if you try to treat it as a physical wavefunction. Neither the energy of this state nor the probability to measure a particular energy eigenstate are well-defined.
Also notice that even if these did converge, they have the wrong dimension! Both expressions have dimensions of energy over length (well your first expression has $1/L^2$, but that should be $1/L$), rather than just energy. So even if they converged, it should be obvious they're not going to give you a physically reasonable result.
A: The difficulties here are caused by the problem in normalizing the position eigenstates
$$
\hat{x}\psi(x)=x_0\psi(x)\longrightarrow \psi(x)=\delta(x-x_0)???\text{ or } |\psi(x)|^2=\delta(x-x_0)???
$$
Since these eigenstates are poorly defined, the straighforward math produces nonsencical results.
So I suggest a different approach: taking an integrable wave function. Indeed,  delta-function is just a limiting case of a series of functions with diminishing width. One can choose many different sequences as a representation (Gaussian, Lorentzian, etc.), but we have here an additional constraint in the form of the boundary conditions at the walls of the well. The analytical choices satisfying this condition may be difficult to work with, but one could probably do with something like:
$$
\psi(x)
=
\begin{cases}
\frac{1}{\sqrt{2\Delta}}, |x-x_0|<\Delta,\\
0, |x-x_0|>\Delta
\end{cases}
$$
One can then calculate all the Fourier coefficients and the necessary quantities and in the end take the limit $\Delta\longrightarrow 0$.
