Fix BEM solving tolerance

[larsoner]

1. Our rhoend was too tight -- it's within the numerical noise of the objective function minimum so on latest SciPy with newer PyPRIMA cobyla it just bounces around until maximum function evals (2000 I think?) hits, which takes about a second. This makes it take < 100 ms for the fit.
2. Cache the fit using relative radii. The values only depend on relative radii, so operate in that space and cache the result. (This speeds up tests but could also speed up end-user cases. The results are quite small so no risk of cache blowing up I think.)
3. Adjust one test whose value changes very slightly due to fit change
4. Add comments justifying rhobeg and rhoend limits, as well as solver choice (I had Claude investigate COBYQA and other solvers and rhobeg/rhoend, our choices are good!)

Part of #14063. Investigated speedups and drafted changes using Claude Fable 5 and Opus 4.8, then I revised (and understand + vouch for all changes).