Parameter Study & Benchmarking Results for solve_poisson_bvp and solve_poisson_ivp

@marco-2023 Heres the detailed parameter study and benchmark result that was asked for.

Parameter Study & Benchmarking Results

To establish reliable default parameters and compare the performance of the BVP and IVP Poisson solvers, I benchmarked both solvers against a normalized s-type Gaussian charge density:

$$
\rho(r)=\left(\frac{\alpha}{\pi}\right)^{3/2}e^{-\alpha r^2}
$$

which has the analytical electrostatic potential:

$$
\phi(r)=\frac{\mathrm{erf}(\sqrt{\alpha}r)}{r}
$$

The benchmark script used to generate these results is located on my fork branch --> here

---

1. BVP Solver Benchmarks (solve_poisson_bvp)

Setup

• Radial Grid: GaussLegendre(n)
• Transform: BeckeRTransform(1e-5, R=1.5)
• Geometry: Single center at the origin
• Angular Degree: 29

---

Study 1 — Effect of Radial Grid Size (n_radial)

Parameters

• tol = 1e-6
• remove_large_pts = 10.0
• α = 0.5

n_radial	Relative L2 Error	Max Absolute Error	Time (s)	Status
30	2.54e-04	1.85e-04	4.2	OK
50	3.72e-05	2.69e-05	4.7	OK
80	6.21e-06	4.47e-06	5.0	OK
100	2.72e-06	1.95e-06	5.0	OK
150	5.41e-07	3.88e-07	5.0	OK
200	1.72e-07	1.23e-07	5.4	OK

Observations

• Approximately fourth-order convergence is observed.
• n_radial = 100 provides the best trade off between accuracy and computational cost.
• Increasing the radial grid beyond 100 points yields diminishing returns.

---

Study 2 — Effect of ODE Solver Tolerance (tol)

Parameters

• n = 100
• remove_large_pts = 10.0
• α = 0.5

tol	Relative L2 Error	Max Absolute Error	Time (s)	Status
1e-2	3.37e-06	2.46e-06	4.2	OK
1e-3	3.37e-06	2.46e-06	4.4	OK
1e-4	2.64e-06	1.95e-06	4.5	OK
1e-5	2.73e-06	1.96e-06	4.5	OK
1e-6 (default)	2.72e-06	1.95e-06	4.4	OK
1e-7	2.72e-06	1.95e-06	4.5	OK

Observations

• Tightening the tolerance below 1e-4 has no measurable effect on accuracy.
• The error floor is determined by the radial grid rather than the ODE solver.
• The current default tol = 1e-6 is safe, although 1e-4 would offer a slight speedup with essentially identical accuracy.

---

Study 3 — Effect of remove_large_pts

Parameters

• n = 100
• tol = 1e-6
• α = 0.5

remove_large_pts	Relative L2 Error	Max Absolute Error	Time (s)	Status
2.0	3.37e-02	2.35e-02	4.6	OK
5.0	6.58e-05	1.09e-04	4.7	OK
10.0	2.72e-06	1.95e-06	4.6	OK
50.0	4.01e-06	2.81e-06	4.8	OK
100.0	4.19e-06	2.93e-06	4.7	OK
1e6 (default)	4.41e-06	3.07e-06	4.7	OK

Observations

• Changing the default threshold from 1e6 to 10.0 reduces relative L2 error by approximately 38% (meaning the default 1e6 has ~62% higher error).
• The improvement results from eliminating far-field asymptotic numerical noise.
• No measurable runtime penalty is observed.

---

Study 4 — BVP Robustness vs Density Sharpness (α)

Parameters

• n = 150
• tol = 1e-6
• remove_large_pts = 10.0

α	Description	Relative L2 Error	Max Absolute Error	Status
0.1	Very wide	1.27e-04	1.23e-04	OK
0.5	Moderate	5.41e-07	3.88e-07	OK
1.0	Standard	5.06e-07	5.07e-07	OK
5.0	Sharp	5.47e-07	1.20e-06	OK
10.0	Nuclear-like	6.56e-07	2.02e-06	OK
50.0	Very sharp	1.41e-06	9.69e-06	OK
100.0	Extremely sharp	2.22e-06	2.16e-05	OK
500.0	Cusp-like	7.85e-06	1.71e-04	OK

Observations

• The BVP solver remains stable across all tested density sharpness values.
• Even for cusp like densities (α = 500), the solver maintains sub-millihartree relative L2 error.

---

2. IVP Solver Benchmarks (solve_poisson_ivp)

Setup

• Radial Grid: Trapezoidal(n)
• Transform: LinearFiniteRTransform(1e-3, 1000.0)
• Integration Interval: (1000.0, 1e-3)

---

Study 5 — IVP Grid Size Requirement

Parameters

• α = 0.5

n_radial	Relative L2 Error	Max Absolute Error	Time (s)	Status
500	5.25e+01	5.31e+01	1.6	Diverged
1000	6.34e-01	8.32e-01	1.6	OK
2000	1.83e-02	3.23e-02	2.6	OK
5000	3.84e-03	1.04e-02	4.7	OK
10000	1.31e-03	4.97e-03	12.9	OK

Observations

• Achieving a relative L2 error of approximately 1.31 × 10⁻³ requires 10,000 radial points.
• At comparable computational scale, the BVP solver (at n_radial = 100) is approximately 500× more accurate while running 2.5× faster than the IVP solver (at n_radial = 10000).

---

Study 6 — IVP Robustness vs Density Sharpness (α)

Parameters

• n = 10000

α	IVP Relative L2	BVP Relative L2	IVP/BVP Error Ratio	Status
0.1	2.68e-03	1.27e-04	21×	OK
0.5	1.31e-03	5.41e-07	2,420×	OK
1.0	2.20e-03	5.06e-07	4,350×	OK
5.0	6.48e-03	5.47e-07	11,800×	OK
10.0	1.47e-02	6.56e-07	22,400×	OK
50.0	1.64e-01	1.41e-06	116,000×	OK
100.0	8.23e-01	2.22e-06	371,000×	OK
500.0	1.09e+01	7.85e-06	1,390,000×	Diverged

Observations

• The IVP solver rapidly loses accuracy as the density becomes increasingly localized.
• At α = 50, the relative L2 error reaches approximately 16%.
• At α = 500, the solution has completely diverged (error > 10.0), whereas the BVP solver remains stable with a relative L2 error of 7.85 × 10⁻⁶.

---

Solver Recommendation

• Prefer solve_poisson_bvp since Its more accurate, more robust, and significantly more computationally efficient than the IVP.
• It integrates naturally with the existing BeckeRTransform grids.

IVP Solver

• It requires extremely dense grids (10,000+ radial points) and becomes unreliable for sharply localized densities.

Suggested Default Update

Update the default value in solve_poisson_bvp:

remove_large_pts = 10.0

Please let me know whether this study and the proposed default parameter updates align with your expectations.