Skip to content

Commit

Permalink
DOC: explain integral examples
Browse files Browse the repository at this point in the history
  • Loading branch information
@redeboer
redeboer committed Feb 11, 2024
1 parent 1239fe7 commit df08529
Showing 1 changed file with 86 additions and 161 deletions.
247 changes: 86 additions & 161 deletions docs/usage/sympy.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,14 +217,37 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Integral"
"## Numerical integrals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Integral over $x^p$ with symbolic range from a to b\n"
"In hadron physics and high-energy physics, it often happens that models contain integrals that do not have an analytical solution.. They can arise in theoretical models, complex scattering problems, or in the analysis of experimental data. In such cases, we need to resort to numerical integrations.\n",
"\n",
"SymPy provides the [`sympy.Integral`](https://docs.sympy.org/latest/modules/integrals/integrals.html#sympy.integrals.integrals.Integral) class, but this does not give us control over whether or not we want to avoid integrating the class analytically. An example of such an analytically unsolvable integral is shown below. Note that the integral does not evaluate despite the `doit()` call."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sp\n",
"\n",
"x, a, b = sp.symbols(\"x a b\")\n",
"p = sp.Symbol(\"p\", positive=True)\n",
"integral_expr = sp.Integral(sp.exp(x) / (x**p + 1), (x, a, b))\n",
"integral_expr.doit()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For amplitude models that contain such integrals that should not be solved analytically, AmpForm provides the {class}`.UnevaluatableIntegral` class. It functions in the same way as [`sympy.Integral`](https://docs.sympy.org/latest/modules/integrals/integrals.html#sympy.integrals.integrals.Integral), but prevents the class from evaluating at all, even if the integral can be solved analytically."
]
},
{
Expand All @@ -239,13 +262,27 @@
},
"outputs": [],
"source": [
"import sympy as sp\n",
"\n",
"from ampform.sympy import UnevaluatableIntegral\n",
"\n",
"x, p, a, b = sp.symbols(\"x,p,a,b\")\n",
"integral_expr = UnevaluatableIntegral(x**p, (x, a, b))\n",
"integral_expr.doit()"
"UnevaluatableIntegral(x**p, (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(x**p, (x, a, b)).doit()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This allows {class}`.UnevaluatableIntegral` to serve as a placeholder in expression trees that we call `doit` on when lambdifying to a numerical function. The resulting numerical function takes **complex-valued** and **multidimensional arrays** as function arguments.\n",
"\n",
"In the following, we see an example where the parameter $p$ inside the integral gets an array as input."
]
},
{
Expand All @@ -260,6 +297,7 @@
},
"outputs": [],
"source": [
"integral_expr = UnevaluatableIntegral(sp.exp(x) / (x**p + 1), (x, a, b))\n",
"integral_func = sp.lambdify(args=[p, a, b], expr=integral_expr)"
]
},
Expand All @@ -277,11 +315,11 @@
"source": [
"import numpy as np\n",
"\n",
"value_a = 1.2\n",
"value_b = 3.6\n",
"value_p = np.array([0.4, 0.6, 0.8])\n",
"a_val = 1.2\n",
"b_val = 3.6\n",
"p_array = np.array([0.4, 0.6, 0.8])\n",
"\n",
"areas = integral_func(value_p, value_a, value_b)\n",
"areas = integral_func(p_array, a_val, b_val)\n",
"areas"
]
},
Expand All @@ -306,186 +344,73 @@
"%config InlineBackend.figure_formats = ['svg']\n",
"import matplotlib.pyplot as plt\n",
"\n",
"x_area = np.linspace(value_a, value_b, num=100)\n",
"x_area = np.linspace(a_val, b_val, num=100)\n",
"x_line = np.linspace(0, 4, num=100)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$x^p$\")\n",
"ax.plot(x_line, x_line ** value_p[0], label=\"$p=0.4$\", color=\"C0\")\n",
"ax.fill_between(x_area, x_area ** value_p[0], alpha=1, color=\"C0\")\n",
"\n",
"ax.plot(x_line, x_line ** value_p[1], label=\"$p=0.6$\", color=\"C1\")\n",
"ax.fill_between(x_area, x_area ** value_p[1], alpha=0.6, color=\"C1\")\n",
"\n",
"ax.plot(x_line, x_line ** value_p[2], label=\"$p=0.8$\", color=\"C2\")\n",
"ax.fill_between(x_area, x_area ** value_p[2], alpha=0.7, color=\"C2\")\n",
"for i, p_val in enumerate(p_array):\n",
" ax.plot(x_line, x_line**p_val, label=f\"$p={p_val}$\", c=f\"C{i}\")\n",
" ax.fill_between(x_area, x_area**p_val, alpha=(0.7 - i * 0.2), color=\"C0\")\n",
"\n",
"ax.text(\n",
" (value_a + value_b) / 2,\n",
" ((value_a ** value_p[0] + value_b ** value_p[0]) / 2) * 0.5,\n",
" \"Area\",\n",
" x=(a_val + b_val) / 2,\n",
" y=((a_val ** p_array[0] + b_val ** p_array[0]) / 2) * 0.5,\n",
" s=\"Area\",\n",
" horizontalalignment=\"center\",\n",
" verticalalignment=\"center\",\n",
")\n",
"ax.annotate(\n",
" \"a\", (value_a, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
")\n",
"ax.annotate(\n",
" \"b\", (value_b, 0.08), textcoords=\"offset points\", xytext=(0, -15), ha=\"center\"\n",
")\n",
"text_kwargs = dict(ha=\"center\", textcoords=\"offset points\", xytext=(0, -15))\n",
"ax.annotate(\"a\", (a_val, 0.08), **text_kwargs)\n",
"ax.annotate(\"b\", (b_val, 0.08), **text_kwargs)\n",
"\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are numerous functions whose integrals cannot be expressed in terms of known special functions and for which no closed-form solutions have been identified yet. Here are a few examples:\n",
"\n",
"1. **Integral Involving Complex Exponential and Polynomials:**\n",
" $$ \\int e^{x + e^{-x}} dx $$\n",
" This integral combines an exponential function with a nested exponential term, leading to a complexity that resists simplification in terms of known functions.\n",
"\n",
"2. **Integral of a Rational Function Times an Exponential:**\n",
" $$ \\int \\frac{e^x}{x^2 + 1} dx $$\n",
" This integral involves the exponential function and a non-trivial rational function, which together do not yield to standard integration techniques.\n",
"\n",
"3. **Integral of an Exponential of a Trigonometric Function:**\n",
" $$ \\int e^{\\sin x} dx $$\n",
" The combination of an exponential function and a sinusoidal function is another case where known methods of integration in terms of elementary or special functions fail.\n",
"\n",
"4. **Product of a Polynomial and a Trigonometric Function:**\n",
" $$ \\int x^4 \\sin(x^2) dx $$\n",
" This integral combines a polynomial with a trigonometric function in a non-standard way, making it resistant to analytical methods.\n",
"\n",
"5. **Integral of a Logarithmic Function Inside an Exponential:**\n",
" $$ \\int e^{\\ln(x)/x} dx $$\n",
" The combination of a logarithmic function inside an exponential function, as seen here, creates a complex relationship not easily resolved by known functions.\n",
"\n",
"In the context of hadron physics and high-energy physics, encountering such challenging integrals is not uncommon. They can arise in theoretical models, complex scattering problems, or in the analysis of experimental data. In such cases, physicists often resort to numerical integration, computational approximations, or the development of new special functions or methods tailored to the problem at hand. These integrals, while not solvable in a traditional sense, can still provide significant insights into the behaviors and properties of the systems being studied."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some example test to compare ``sp.Integral`` and ``UnevaluatableIntegral``"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"UnevaluatableIntegral(sp.exp(-(x**4)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(sp.exp(-(x**4)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"UnevaluatableIntegral(sp.sin(x) / x, (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(sp.sin(x) / x, (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"UnevaluatableIntegral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(sp.exp(x + sp.exp(-x)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"UnevaluatableIntegral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(sp.exp(sp.log(x) / x), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"UnevaluatableIntegral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sp.Integral(sp.exp(sp.sin(x)), (x, a, b)).doit()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"UnevaluatableIntegral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
"The arrays can be complex-valued as well. This is particularly useful when calculating dispersion integrals (see **[TR-003](https://compwa.github.io/report/003#general-dispersion-integral)**)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"sp.Integral(sp.exp(x) / (x**2 + 1), (x, a, b)).doit()"
"integral_func(\n",
" p=np.array([1.5 - 8.6j, -4.6 + 5.5j]),\n",
" a=a_val,\n",
" b=b_val,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## Summations"
]
Expand Down

0 comments on commit df08529

Please sign in to comment.