{
"cells": [
{
"cell_type": "markdown",
"id": "53317fc2",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# A.6 The `scipy` package for scientific computing"
]
},
{
"cell_type": "markdown",
"id": "504e7653",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
" The three main Python packages for exploratory computing are `numpy`, `matplotlib`, and `scipy`. The `numpy` and `matplotlib` packages were already discussed in previous examples. Here, some of the functionalities of `scipy` are demonstrated, including root-finding, minimization, special functions, and numerical integration."
]
},
{
"cell_type": "markdown",
"id": "733f9ced",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## SciPy optimize functions\n",
"\n",
" The function `fsolve` from the `scipy.optimize` package can be used to find the root or roots of a function. In the code cell below, the root of the function $f(x) = \\text{e}^{x/5} - 2$ is found numerically using the `fsolve` function (it can, for this case, also be computed analytically, of course). A starting value for the numerical search must be supplied. The closer the starting value to the root, the more likely (and the quicker) `fsolve` finds the root."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "fda5c5f5",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from myst_nb import glue\n",
"%config InlineBackend.rc = {'figure.dpi': 300}\n",
"%config InlineBackend.figure_format = 'svg'"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "115e8510",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"root=[3.4657359], function value=[2.7533531e-14]\n"
]
}
],
"source": [
"def ffunc(x):\n",
" return np.exp(x / 5) - 2\n",
"\n",
"from scipy.optimize import fsolve\n",
"xroot = fsolve(ffunc, x0=0) # start searching at x0=0\n",
"print(f'root={xroot}, function value={ffunc(xroot)}')"
]
},
{
"cell_type": "markdown",
"id": "ee7be8c3",
"metadata": {},
"source": [
" Next, three points are found where $f(x)$ intersects the function $g(x) = \\sin(x)$. The function `fdiff` computes the difference between $f(x)$ and $g(x)$, and `fsolve` is used to find three roots of `fdiff` for three different starting locations; that is, `fsolve` finds three points where the two functions intersect."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "9980b662",
"metadata": {
"tags": [
"remove-output"
]
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def gfunc(x):\n",
" return 3 * np.sin(x)\n",
"\n",
"def fdiff(x):\n",
" return ffunc(x) - gfunc(x)\n",
"\n",
"x_intersect = fsolve(fdiff, np.array([3, 6, 9])) # 3 starting points\n",
"x = np.linspace(0, 3 * np.pi)\n",
"plt.plot(x, ffunc(x), label=\"f(x)\")\n",
"plt.plot(x, gfunc(x), label=\"g(x)\")\n",
"plt.plot(x_intersect, ffunc(x_intersect), 'k.')\n",
"plt.legend();"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "cac9d8cc",
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [
{
"data": {
"application/papermill.record/image/svg+xml": "\n\n\n",
"application/papermill.record/text/plain": ""
},
"metadata": {
"scrapbook": {
"mime_prefix": "application/papermill.record/",
"name": "intersect-plot"
}
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(4, 2.5))\n",
"plt.plot(x, ffunc(x), label=\"f(x)\")\n",
"plt.plot(x, gfunc(x), label=\"g(x)\")\n",
"plt.plot(x_intersect, ffunc(x_intersect), 'k.')\n",
"plt.legend()\n",
"glue(\"intersect-plot\", fig, display=False)"
]
},
{
"cell_type": "markdown",
"id": "c581f7df",
"metadata": {},
"source": [
"```{glue:figure} intersect-plot\n",
":name: fig_intersect_plot_A.6\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "db7c0f8a",
"metadata": {},
"source": [
" The minimum of a function may be found with the `fmin` function of `scipy.optimize`. The code cell below finds the $x$-location where the difference between the functions $g(x)$ and $h(x) = \\text{e}^{x/5} + 2$ is the smallest, starting the search at $x = 1.5$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9689b132",
"metadata": {
"tags": [
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.356826\n",
" Iterations: 11\n",
" Function evaluations: 22\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def hfunc(x):\n",
" return np.exp(x/5) + 2\n",
"\n",
"def fdiff(x):\n",
" return hfunc(x) - gfunc(x)\n",
"\n",
"from scipy.optimize import fmin\n",
"xmin = fmin(fdiff, 1.5)\n",
"\n",
"plt.plot(x, hfunc(x), label=\"h(x)\")\n",
"plt.plot(x, gfunc(x), label=\"g(x)\")\n",
"plt.plot(xmin, hfunc(xmin), 'C0.')\n",
"plt.plot(xmin, gfunc(xmin), 'C1.')\n",
"plt.legend();"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "2feec6f9",
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [
{
"data": {
"application/papermill.record/image/svg+xml": "\n\n\n",
"application/papermill.record/text/plain": ""
},
"metadata": {
"scrapbook": {
"mime_prefix": "application/papermill.record/",
"name": "min-plot"
}
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(4, 2.5))\n",
"plt.plot(x, hfunc(x), label=\"h(x)\")\n",
"plt.plot(x, gfunc(x), label=\"g(x)\")\n",
"plt.plot(xmin, hfunc(xmin), 'C0.')\n",
"plt.plot(xmin, gfunc(xmin), 'C1.')\n",
"plt.legend()\n",
"glue(\"min-plot\", fig, display=False)"
]
},
{
"cell_type": "markdown",
"id": "c0233235",
"metadata": {},
"source": [
"```{glue:figure} min-plot\n",
":name: fig_min_plot_A.6\n",
"```"
]
},
{
"cell_type": "markdown",
"id": "99a23fa8",
"metadata": {},
"source": [
"## SciPy special functions\n",
" Two other subpackages of the `scipy` package are `scipy.special` for special functions and `scipy.integrate` for numerical integration. As an example, consider the exponential integral defined as\n",
"\n",
"$$\n",
"\\text{E}_1(u) = \\int_u^{\\infty} \\frac{e^{-s}}{s} ds \n",
"$$\n",
"\n",
"$\\text{E}_1$ is available as `exp1` from the `scipy.special` package, but can also be evaluated using numerical integration. The `quad` function of the `scipy.integrate` package is used and takes as input the function to integrate and the lower and upper limits of the integration; use `np.inf` for $\\infty$. Note that `quad` returns two values: the outcome of the numerical integration and an error estimate. The results of the numerical integration and the special function are compared below for $\\text{E}_1(1)$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4a4fa1c9",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E1(1) from numerical integration: 0.2193839343955124\n",
"E1(1) from special function exp1: 0.2193839343955205\n"
]
}
],
"source": [
"from scipy.integrate import quad\n",
"from scipy.special import exp1\n",
"\n",
"def integrand(s):\n",
" return np.exp(-s) / s\n",
"\n",
"E1num = quad(integrand, 1, np.inf)[0] # numerical integration\n",
"\n",
"print('E1(1) from numerical integration:', E1num)\n",
"print('E1(1) from special function exp1:', exp1(1))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}