  Boundary element methods produce dense linear systems that can be accelerated
via multipole expansions. Solved with Krylov methods, this implies computing
the matrix-vector products within each iteration with some error, at an
accuracy controlled by the order of the expansion, $p$. We take advantage of a
unique property of Krylov iterations that allow lower accuracy of the
matrix-vector products as convergence proceeds, and propose a relaxation
strategy based on progressively decreasing $p$. In extensive numerical tests of
the relaxed Krylov iterations, we obtained speed-ups of between $2.1\times$ and
$3.3\times$ for Laplace problems and between $1.7\times$ and $4.0\times$ for
Stokes problems. We include an application to Stokes flow around red blood
cells, computing with up to 64 cells and problem size up to 131k boundary
elements and nearly 400k unknowns. The study was done with an in-house
multi-threaded C++ code, on a hexa-core CPU. The code is available on its
version-control repository,
\href{https://github.com/barbagroup/fmm-bem-relaxed}{https://github.com/barbagroup/fmm-bem-relaxed}.
