DOC: explain implementationm

Author: redeboer

.cspell.json

@@ -266,6 +266,7 @@
     "isospin",
     "itertools",
     "jupyter",
+    "Kronrod",
     "Källén",
     "lambdification",
     "lambdified",

docs/usage/dynamics/integration-algorithms.ipynb

@@ -78,6 +78,27 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "The dispersion integral of the phasespace factor cannot be solved analytically for higher angular momenta ($\\ell > 0$). We therefore need to perform a numerical integration. However, since the integrand of the dispersion integral contains singularities, it is important to choose the numerical integration algorithm.\n",
+    "\n",
+    "This notebook contains an interactive widget that shows the numerical stability and performance of different numerical integration algorithms when computing the dispersion integral along (close to) the real axis of the complex energy plane.\n",
+    "\n",
+    ":::{tip}\n",
+    "The conclusion is that [Romberg's method](https://en.wikipedia.org/wiki/Romberg%27s_method) ({func}`quadax.romberg`) is the best choice for computing the dispersion integral close to the physical (real) axis. For computing the dispersion integral further away from the axis, the [Gauss-Kronrod method](https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula) (e.g. {func}`quadax.quadgk`) are accurate enough and faster than Romberg's method.\n",
+    "\n",
+    "- Use {func}`quadax.romberg` when computing the dispersion integral for an amplitude model fit (physical axis).\n",
+    "- Use {func}`quadax.quadgk` when determining the pole positions of the amplitude model (complex plane).\n",
+    "\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Chew–Mandelstam dispersion integral\n",
+    "\n",
+    "The general form of the Chew–Mandelstam dispersion integral $\\Sigma_\\ell(s)$ is an integral over the product $\\rho(s) \\, n_\\ell^2(s)$ of the phase space factor $\\rho(s)$ ({obj}`.PhaseSpaceFactor`) and the square of the form factor $n_\\ell^2(s)$ ({obj}`.FormFactor`). The required formulas are:\n",
+    "\n",
     "```{autolink-skip}\n",
     "\n",
     "```"
@@ -121,6 +142,43 @@
     "Math(src)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here, $i\\epsilon$ indicates that $\\Sigma_\\ell(s)$ with $s\\in\\mathbb{R}$ is formulated just above the real axis in order to avoid $s'-s=0$. However, as can be seen in the widget below, $i\\epsilon$ is not required when $s\\in\\mathbb{C}$, giving us the general form:\n",
+    "\n",
+    "$$\n",
+    "\\begin{aligned}\n",
+    "\\Sigma_\\ell\\left(s\\right) \\;&=\\; \\frac{s - s_\\mathrm{thr}}{\\pi}\n",
+    "  \\int\\limits_{s_\\mathrm{thr}}^{\\infty}\n",
+    "  \\frac\n",
+    "  {\\rho\\!\\left(s'\\right) n_\\ell^2\\!\\left(s'\\right) ds'}\n",
+    "  {\\left(s' - s_\\mathrm{thr}\\right) \\left(s'- s\\right)} \\\\\n",
+    "n_\\ell^2(s') \\;&=\\; \\mathcal{F}_\\ell^2\\!\\left(s', m_1, m_2\\right) \\\\\n",
+    "s_\\mathrm{thr} \\;&=\\; (m_1 + m_2)^2\n",
+    "\\end{aligned}\n",
+    "$$\n",
+    "\n",
+    "The reason is that the function $\\rho(s')\\,n_\\ell^2(s')$ does not have a branch cut along the real axis above the $s_\\mathrm{thr}$ threshold, so there is no need to approach some contour around such a potential discontinuity. The form of the dispersion integral is explained and derived [here](https://redeboer.github.io/phd-thesis/chapter2#sec-analytic-continuation).\n",
+    "\n",
+    "When using this form of the dispersion integral to compute $\\Sigma_\\ell(s)$ for $s\\in\\mathbb{R}$, the caller should use the fact that $\\lim_{\\epsilon\\to 0^+} \\Sigma_\\ell(s + i\\epsilon)$. In the numerical implementation, that means giving `s + epsilon * 1j` as input when `s` is a array with `float`s. The challenge is to find a value for $\\epsilon$: the smaller the value for $\\epsilon$, the closer we are to the physical axis, but the less accurate the numerical integration becomes. The widget below allows you to explore this trade-off for different numerical integration algorithms."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Numerical implementation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this notebook, we implement the formulas [listed above](#chewmandelstam-dispersion-integral) numerical rather than lambdifying the symbolic expressions, so that we have full control over the numerical implementation and to make the implementation recognizable to pure array-oriented workflows."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -206,9 +264,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the case of S-waves, we can compare the result of the integration to the analytical result.\n",
-    "\n",
-    "Note that in this case, we define the break-up momentum in the same way as in {class}`.BreakupMomentum`, because in the widget below, we only evaluate this function along the real axis, where the cut structure doesn't matter and we prefer performance (see {ref}`usage/dynamics/analytic-continuation:Numerical precision and performance`)."
+    "In the case of $S$-waves ($\\ell=0$), we can compare the result of the integration to the analytical solution to the integral:"
    ]
   },
   {
@@ -253,6 +309,13 @@
     "    )"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interactive visualization"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,