I found a text of Robert Israel about Richardson extrapolation
http://www.math.ubc.ca/~israel/m215/rich/rich.html
There is the formula
$$Q_R=\frac{h_2^pQ(h_1)-h_1^pQ(h_2)}{h_2^p-h_1^p}+\mathcal O(h^{p+1})$$
where $Q=Q(h)+\mathcal O(h^p)$. I wonder if this can be used together with an algorithm that evaluate $a_n$ in some convergent real sequence, letting $h_i=1/n_i$. Normally I don't know $p$ and I wonder if that problem generally can be solved by a trial-and-error method?
Concrete example. If $(s_n)_n$ is the ordered sequence of square free numbers, it can be proved that $a_n=s_n/n$ converges to $\pi^2/6$. Therefor any strictly increasing sequence $(t_n)_n$ of natural numbers that includes all square free numbers have $\displaystyle 1\leq\lim_{n\to\infty}t_n/n\leq\pi^2/6$.
I have chosen to examine the sequence that include all square free numbers and numbers with at most one square prime in the prime factorization. My algorithm works brutal force by testing the numbers $1<k\leq n$ and I don't know the number $p$ in the formula but I first tried $p=1$.
$\displaystyle\frac{t_{1,000,000}}{1,000,000}\approx 1.2365450\;$ and $\displaystyle\frac{t_{10,000,000}}{10,000,000}\approx 1.2365731\;$ and $\displaystyle\frac{t_{100,000,000}}{100,000,000}\approx 1.2365767$.
The formula gives for $p=1$
$\displaystyle\frac{10^{-7}\cdot1.2365450-10^{-6}\cdot1.2365731}{10^{-7}-10^{-6}}\approx 1.2365762\;$ and
$\displaystyle\frac{10^{-8}\cdot1.2365731-10^{-7}\cdot1.2365767}{10^{-8}-10^{-7}}\approx 1.2365771$
And for $p=2$
$\displaystyle\frac{10^{-14}\cdot1.2365450-10^{-12}\cdot1.2365731}{10^{-14}-10^{-12}}\approx 1.2365734\;$ and
$\displaystyle\frac{10^{-16}\cdot1.2365731-10^{-14}\cdot1.2365767}{10^{-16}-10^{-14}}\approx 1.2365767$
Can I from this draw the conclusion that $p=1$ and that the limit is close to $1.2365771$?
No, you cannot discern the order of the principal error term from Richardson's error estimate. You must either derive an asymptotic error expansion or observe the behavior of Richardson's fraction, see below.
Suppose $a_n \rightarrow a$ as $n = \rightarrow \infty$ and $n \in \mathbb{N}$. Then you can apply Richardson's techniques provided the error, i.e., $e_n = a - a_n$, satisfies an asymptotic error expansion of the form $$ e_n = \alpha n^{-p} + \beta n^{-q} + O(n^{-r}), $$ where $\alpha$, $\beta$, $p$, $q$, and $r$ are constants independent of $n$ and $0 < p < q < r$. The exponents $p$, $q$ and $r$ need not be integers. The term $\alpha n^{-p}$ is the primary error term and $p$ is its order.
In this context, it is straightforward to verify that Richardson's fraction $$ F_{n} = \frac{a_{n/2}-a_{n/4}}{a_{n}-a_{n/2}}, \quad n = 2^k$$ converges to $2^p$ as $k \rightarrow \infty$ and $k \in \mathbb{N}$. Moreover, for sufficiently large $k$, the convergence is strictly increasing when $\beta/\alpha < 0$ and strictly decreasing when $\beta/\alpha > 0$. Deviation from this pattern is due to the fact that the higher order terms are still significant (real and machine arithmetic) or due to the fact that subtractive cancellation can no longer be ignored when computing the nominator and denominators for $F_n$ (machine arithmetic only). Within the asymptotic range where the computed value of Richardson's fraction behaves as the exact value you can really on Richardson's error estimate, i.e., $$ a - a_n = \frac{a_{n} - a_{n/2}}{2^p -1} + O(n^{-q}).$$ The accuracy of Richardson's error estimate tend to increase until subtractive cancellation sets in.
The Wikipedia page on square-free integers strongly suggest that the order of the primary error term $p$ is not an integer.