Let's consider this boundary-value problem:
$\begin{cases} -\Delta V = \rho & \rm{in}\ \Omega \\ V=0 & \rm{on}\ \partial \Omega \end{cases}.$
We know that this problem has a variational formulation: its solutions are critical points of the energy functional
$I[u]=\int_{\Omega} \left( \frac{1}{2}\lvert \nabla u \rvert^2 - u\rho \right)\, dx,$
(cfr. Evans Partial differential equations 2nd ed., §2.2.5).
Now let us look at one (of the many) possible physical interpretations of this problem, the electrostatic one. If units are chosen so that $\varepsilon_0=1$, this equation describes the electric potential $V$ in a region $\Omega$ filled with a charge distribution $\rho$ and whose boundary is grounded.
It would be nice if the energy functional introduced above coincided with the total energy of this physical system. Unfortunately, it seems to me not to be so. In fact, as I read in Feynman's Lectures on physics, vol.II, §8-5, the total energy of this system consists of two additive parts: one is
$\frac{1}{2}\int_{\Omega} \lvert \nabla V \rvert^2\, dx$
and is due to the electric field; the other is
$\frac{1}{2}\int_{\Omega} V\rho\, dx$
and is due to the charge distribution.
This induces me to think that the energy functional should be
$U[u]=\int_{\Omega} \left(\frac{1}{2}\lvert \nabla u \rvert^2 + \frac{1}{2}u\rho \right)\, dx,$ instead of the correct $I[u]=\int_{\Omega} \left( \frac{1}{2}\lvert \nabla u \rvert^2 - u\rho \right)\, dx.$ Why am I wrong?