  This paper analyzes the influence of general, small volume, inclusions on the
trace at the domain's boundary of the solution to elliptic equations of the
form $\nabla \cdot D^\eps \nabla u^\eps=0$ or $(-\Delta + q^\eps) u^\eps=0$
with prescribed Neumann conditions. The theory is well-known when the
constitutive parameters in the elliptic equation assume the values of different
and smooth functions in the background and inside the inclusions. We generalize
the results to the case of arbitrary, and thus possibly rapid, fluctuations of
the parameters inside the inclusion and obtain expansions of the trace of the
solution at the domain's boundary up to an order $\eps^{2d}$, where $d$ is
dimension and $\eps$ is the diameter of the inclusion. We construct inclusions
whose leading influence is of order at most $\eps^{d+1}$ rather than the
expected $\eps^d$. We also compare the expansions for the diffusion and
Helmholtz equation and their relationship via the classical Liouville change of
variables.
