We seek approximations of the solution u of the Neumann problem for the equation
L
u
=
f
Lu = f
in
Ω
\Omega
with special emphasis on high-order accuracy at a given point
x
0
∈
Ω
¯
{x_0} \in \bar \Omega
. Here
Ω
\Omega
is a bounded domain in
R
N
(
N
⩾
2
)
{R^N}(N \geqslant 2)
with smooth boundary, and L is a second-order, uniformly elliptic, differential operator with smooth coefficients. An approximate solution
u
h
{u_h}
is determined by the standard Galerkin method in a space of continuous piecewise polynomials of degree at most
r
−
1
r - 1
on a partition
Δ
h
(
x
0
,
α
)
{\Delta _h}({x_0},\alpha )
of
Ω
\Omega
. Here h is a global mesh-size parameter, and
α
\alpha
is the degree of a certain systematic refinement of the mesh around the given point
x
0
{x_0}
, where larger
α
\alpha
’s mean finer mesh, and
α
=
0
\alpha = 0
corresponds to the quasi-uniform case with no refinement. It is proved that, for suitable (sufficiently large)
α
\alpha
’s the high-order error estimate
(
u
−
u
h
)
(
x
0
)
=
O
(
h
2
r
−
2
)
(u - {u_h})({x_0}) = O({h^{2r - 2}})
holds. A corresponding estimate with the same order of convergence is obtained for the first-order derivatives of
u
−
u
h
u - {u_h}
. These estimates are sharp in the sense that the required degree of refinement in each case is essentially the same as is needed for the local approximation to this order near
x
0
{x_0}
. For the estimates to hold, it is sufficient that the exact solution u have derivatives to the rth order which are bounded close to
x
0
{x_0}
and square integrable in the rest of
Ω
\Omega
. The proof of this uses high-order negative-norm estimates of
u
−
u
h
u - {u_h}
. The number of elements in the considered partitions is of the same order as in the corresponding quasi-uniform ones. Applications of the results to other types of boundary value problems are indicated.