{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerical Integration\n", "\n", "## Overview, Objectives, and Key Terms\n", " \n", "In [Numerical Differentiation](L19_NumericalDifferentiation.ipynb), finite-difference approximations were developed and used to evaluate derivatives numerically. In this lecture, we'll explore and implement our own numerical methods for evaluation of integrals. We'll also continue using SciPy, this time through its numerical integration tools.\n", " \n", "### Objectives\n", "\n", "By the end of this lesson, you should be able to\n", "\n", "- Evaluate definite integrals numerically using left- and right-sided Riemann sums, the mid-point rule, and the trapezoid rule.\n", "- Evaluate integrals of sampled data using `scipy.integrate.trapezoid`, `scipy.integrate.simpson`, and `scipy.integrate.cumulative_trapezoid`.\n", "- Interpret how sampling choices create overestimates and missed area.\n", "- Evaluate integrals of callable functions using `scipy.integrate.quad` and related functions.\n", "- Establish the order of an integration scheme using numerical, graphical, or symbolic means.\n", "- Recognize how numerical integration can be used to build probability and balance checks in deterministic response models.\n", "\n", "### Key Terms\n", "\n", "- Riemann sum\n", "- left-sided Riemann sum\n", "- right-sided Riemann sum\n", "- mid-point rule\n", "- trapezoid rule\n", "- Simpson's rule\n", "- sampled-data integration\n", "- adaptive quadrature\n", "- `scipy.integrate.trapezoid`\n", "- `scipy.integrate.simpson`\n", "- `scipy.integrate.cumulative_trapezoid`\n", "- `scipy.integrate.quad`\n", "- `scipy.integrate.dblquad`\n", "- $n$th-order approximation\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Riemann Sums in Practice\n", "\n", "For derivatives, we took a limit of a finite difference for an exact answer but stopped that limit for a sufficiently small $\\Delta$ to ensure a reasonable estimate. For the first derivative, we saw the first-order forward and backward schemes along with the second-order central difference scheme. \n", "\n", "For integrals, we have analogs to these choices. You may remember the following formal definition for a definite integral:\n", "\n", "$$\n", " \\int^b_a f(x) dx = \n", " \\lim_{n \\to \\infty} \\sum^{n}_{i=1} f(i (b-a)/n + a) (b-a)/n\n", "$$\n", "\n", "Here, the sum is a *Riemann sum* (specifically a right-sided sum). Integrals are areas under curves, and Riemann sums approximate those areas by rectangles. Just as for derivatives, we *don't take the limit* but let $n$ be a finite number. Moreover, we have different numerical options depending on how we define the rectangles of the Riemann sum. \n", "\n", "Let's look at this concretely for $f(x) = x^2$ over $x\\in[0, 1]$ with $n=5$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "For this increasing function, the right-sided estimate has positive signed error: the estimate lies above the curve on each interval." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the height of a rectangle is equal to the function evaluated at the rectangle's right $x$ coordinate. We call such a sum a *right-sided* Riemann sum. For this particular function, this right-sided sum leads to some significant over estimation. \n", "\n", "If instead we defined our sum to be\n", "\n", "$$\n", " S =\\frac{b-a}{n} \\sum^{n}_{i=1} f \\Bigg ( \\frac{ (i-1)(b-a) }{n} + a \\Bigg ) \n", "$$\n", "\n", "the picture becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The left-sided estimate has negative signed error here: the estimate misses area under the curve on each interval." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the rectangle height is based on the function evaluated at its \n", "left side, making this a *left-sided* Riemann sum. Just as the right-sided\n", "sum *overestimates* the integral, the left-sided sum *underestimates* \n", "the integral." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Exercise**: Given the points $x = [0, 2/3, 4/3, 2]$ and corresponding values $y = [1, 3/5, 3/7, 1/3]$, estimate $\\int^2_0 y(x)dx$ using the right-sided sum and left-sided sum. What is the true value (assuming $y = (x+1)^{-1}$)?\n", "\n", "
\n", "\n", "**Exercise**: Consider the following two columns of data:\n", "\n", "| $x$ | $\\sqrt{x}$ |\n", "|-------|------------|\n", "|0.0000 | 0.0000 |\n", "|0.3333 | 0.5774 |\n", "|1.0000 | 1.0000 |\n", "|1.5000 | 1.2247 |\n", "|2.0000 | 1.4142 |\n", "\n", "Estimate $\\int^2_0 f(x) dx$ using the right-sided sum and left-sided sum. Note that the $x$ points are not evenly spaced.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Beyond the One-Side Sums\n", "\n", "Can we do better? Remember, I said that the forward, backward, and central\n", "finite-difference approximations have analogous Riemann sums. Here,\n", "take the right- and left-sided sums to be like the forward and backward\n", "differences. If you hadn't noticed, the central difference is nothing\n", "but an *average* of the forward and backward differences:\n", "\n", "$$\n", " \\frac{1}{2} \\Bigg ( \\frac{f(x+h)-f(x)}{h} + \\frac{f(x)-f(x-h)}{h} \\Bigg ) = \\frac{f(x+h)-f(x-h)}{2h}\n", "$$\n", "\n", "so we might expect something similar for the Riemann sums. There are\n", "two possibilities. First, we can evaluate the rectangle heights at the \n", "average value of its $x$ coordinates, i.e.,\n", "\n", "$$\n", " S = \\frac{b-a}{n} \\sum^{n}_{i=1} f \\Bigg ( \\frac{ (i-1/2)(b-a) }{n} + a \\Bigg ) \n", "$$\n", " \n", "for which the picture becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The midpoint estimate mixes small positive and negative error regions. For this smooth curve, those regions cancel much more effectively." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This method is often called the **midpoint rule**, and it's probably my \n", "favorite integration scheme because it's pretty accurate and dirt simple\n", "to implement.\n", "\n", "The alternative is to average the right- and left-sided sums, which leads\n", "to \n", "\n", "$$\n", "\\begin{split}\n", " S &= \\frac{b-a}{2n} \n", " \\Bigg [ \\sum^{n}_{i=1} f \\Bigg ( \\frac{ (i)(b-a) }{n} + a \\Bigg ) \n", " + \\sum^{n}_{i=1} f \\Bigg ( \\frac{ (i-1)(b-a) }{n} + a \\Bigg ) \\Bigg ] \\\\\n", " &= \\frac{b-a}{2n} \n", " \\Bigg [ f(a) + f(b) + \\sum^{n-1}_{i=1} f \\Bigg ( \\frac{ i(b-a) }{n} + a \\Bigg ) \\Bigg ] \n", "\\end{split}\n", "$$\n", " \n", "If you consider this carefully, the resulting geometric picture turns\n", "rectangles into trapezoids (because the average of the area of \n", "two rectangles of heights $y(x_1)$ and $y(x_2)$ that share \n", "the same $x$ coordinates $x_1$ and\n", "$x_2$ is equal to the area of the trapezoid whose coordinates\n", "include $(x_1, y_1)$ and $(x_2, y_2)$),\n", "and the method is known as the **trapezoid rule**. (It's those pesky\n", "end-point terms $(a)$ and $f(b)$\n", "that make this rule just a hair less convenient to \n", "apply in some cases, so I always stick to the midpoint rule when I can!)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The trapezoid estimate uses a straight-line segment across each interval. The colored regions show where that line differs from the curve." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a function $f(x)$ with sufficiently nice properties (beyond\n", "our scope), it's not too difficult to estimate the error associated\n", "with these integrals as $\\Delta = (b-a)/n$ goes to zero (or \n", "equivalently, $n \\to \\infty$. However, I'll skip giving away the \n", "answers and let you determine these properties by numerical experiments\n", "in the homework!\n", "\n", "Finally, the integration rules just explored represent the simplest \n", "possible [Newton-Cotes](http://mathworld.wolfram.com/Newton-CotesFormulas.html) formulas. The basic idea of Newton-Cotes \n", "integration is that a function $f(x)$ is assumed to be piece-wise\n", "polynomial in finite ranges and evaluated at evenly-spaced points $x_i$. Left-sided sums, right-sided sums, and\n", "and the midpoint rule take \n", "$f(x)$ to be *constant* between adjacent points $x_i$ and\n", "$x_{i+1}$, while the trapezoid rule assumes $f(x)$ is \n", "*linear* between adjacent points. If $f(x)$ is assumed to \n", "be quadratic between $x_i$ and $x_{i+2}$, the \n", "result is [Simpson's rule](http://mathworld.wolfram.com/SimpsonsRule.html) (another favorite of mine). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Library Tools for Sampled Data\n", "\n", "The rules above are simple enough that you should be able to implement them yourself. That is still important: if you cannot write a trapezoid rule in a few lines, then the library function is just a magic word.\n", "\n", "Once the idea is clear, though, use the library routine. For data that have already been sampled at points `x`, the most common tools are\n", "\n", "- `scipy.integrate.trapezoid(y, x)` for the trapezoid rule,\n", "- `scipy.integrate.simpson(y, x=x)` for Simpson's rule, and\n", "- `scipy.integrate.cumulative_trapezoid(y, x, initial=0.0)` when you want the accumulated integral as a function of position or time.\n", "\n", "This is a different situation from integrating a function like `np.exp` over `[0,1]`. Here, the data are already on the page. The job is to integrate those samples without pretending we know more about the function than the data actually tell us." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2026-07-02T15:33:03.850936Z", "iopub.status.busy": "2026-07-02T15:33:03.850841Z", "iopub.status.idle": "2026-07-02T15:33:04.013100Z", "shell.execute_reply": "2026-07-02T15:33:04.012699Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "trapezoid = 1.837937\n", "simpson = 1.901891\n", "reference = 1.885618\n", "[0. 0.09621061 0.6220107 1.17819692 1.83793653]\n" ] } ], "source": [ "import numpy as np\n", "from scipy.integrate import cumulative_trapezoid, simpson, trapezoid\n", "\n", "x = np.array([0.0000, 0.3333, 1.0000, 1.5000, 2.0000])\n", "y = np.sqrt(x)\n", "\n", "I_trap = trapezoid(y, x)\n", "I_simp = simpson(y, x=x)\n", "I_cumulative = cumulative_trapezoid(y, x, initial=0.0)\n", "I_ref = (2/3)*2**(3/2)\n", "\n", "print(f\"trapezoid = {I_trap:.6f}\")\n", "print(f\"simpson = {I_simp:.6f}\")\n", "print(f\"reference = {I_ref:.6f}\")\n", "print(I_cumulative)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The cumulative result is often more useful than the final number. For example, if the data were an acceleration history, the cumulative integral would give velocity history; a second cumulative integral would give displacement history. In NE language, the same habit shows up whenever a local quantity has to be accumulated into a total response." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "**Exercise**: Suppose you accessed your fancy phone's accelerometer log and extracted the following accelerations (in m/s$^2$) and times (in s):\n", "```\n", "a = 0, 0.1, 0.3, -0.2, 0.1\n", "t = 2, 4, 6, 8, 10\n", "```\n", "Estimate how far you went in those 8 seconds using the trapezoid rule.\n", "
\n", "\n", "**Exercise**: Recall that the incremental work $\\delta W$ done by an expanding gas is given by $\\delta W = p \\delta V$, where $p$ is pressure and $\\delta V$ is the change in volume. Consider the following table of measured pressures (MPa) for several volumes (cm ^3^) during the expansion of the gas within an internal combustion engine.\n", "\n", "| $p$ | $V$ |\n", "|---------|----------|\n", "| 15.6835 | 79.8508 |\n", "| 12.7127 | 88.6913 |\n", "| 9.3703 | 103.3057 |\n", "| 6.5612 | 123.4549 |\n", "| 4.5031 | 149.0202 |\n", "| 2.0138 | 222.8411 |\n", "\n", "Estimate the work performed by the gas using the trapezoid rule. \n", "\n", "
\n", "\n", "**Exercise**: Revisit the last exercise by determinging the $n$ that satisfies $pV^n = c$ for some constant $c$. Then, using the same measured data, apply the midpoint rule to estimate the work done by the gas.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What About Accuracy?\n", "\n", "As was done for the finite-difference methods, it is possible to determine the *order* of various numerical integration rules. The basic idea is to observe the (absolute) error between an exact, reference value for an integral and the approximation of that integral for various step sizes $\\Delta$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "The black points show the right-sided error for a smooth test function. The purple points show a second, less smooth comparison case; the reference lines give first- and second-order visual guides." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "**Exercise**: Determine whether the right-sided sum is first order or second order using the definite integral $\\int^1_0 e^x dx = e - 1$ as a test case. \n", "\n", "*Solution*:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2026-07-02T15:33:04.016906Z", "iopub.status.busy": "2026-07-02T15:33:04.016813Z", "iopub.status.idle": "2026-07-02T15:33:04.152650Z", "shell.execute_reply": "2026-07-02T15:33:04.151831Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "<>:24: SyntaxWarning: invalid escape sequence '\\m'\n", "<>:25: SyntaxWarning: invalid escape sequence '\\m'\n", "<>:26: SyntaxWarning: invalid escape sequence '\\D'\n", "<>:24: SyntaxWarning: invalid escape sequence '\\m'\n", "<>:25: SyntaxWarning: invalid escape sequence '\\m'\n", "<>:26: SyntaxWarning: invalid escape sequence '\\D'\n", "/tmp/ipykernel_143228/2207420918.py:24: SyntaxWarning: invalid escape sequence '\\m'\n", " plt.loglog(Delta, Delta, 'r--', label='$\\mathcal{O}(\\Delta)$')\n", "/tmp/ipykernel_143228/2207420918.py:25: SyntaxWarning: invalid escape sequence '\\m'\n", " plt.loglog(Delta, 10*Delta**2, 'b-.', label='$\\mathcal{O}(\\Delta^2)$')\n", "/tmp/ipykernel_143228/2207420918.py:26: SyntaxWarning: invalid escape sequence '\\D'\n", " plt.xlabel('$\\Delta$')\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAG1CAYAAAAxwRT8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXBxJREFUeJzt3Qd0lNXWBuA39A7SO6IoRWpoAiJFCEVp0juIKIoFURFEBSwXryJXL00p0ntHQBEBAeFXukoHwQvSQu+hzr/eczKTQkvCTKa9z1qjzJdJ8uVLmT1n77N3iMPhcEBERETETyTx9gmIiIiIxIeCFxEREfErCl5ERETEryh4EREREb+i4EVERET8ioIXERER8SsKXkRERMSvKHgRERERv5IMAebmzZs4fPgw0qdPj5CQEG+fjoiIiMQBe+aeP38euXPnRpIkSYIreGHgki9fPm+fhoiIiCTAwYMHkTdv3uAKXrji4vziM2TI4O3TERERkTg4d+6cWXxwPo8HVfDiTBUxcFHwIiIi4l/iUvKhgl0RERHxKwpeRERExK8EXNoorm7cuIFr1655+zQCRvLkyZE0aVJvn4aIiASBZMG4Fevo0aM4c+aMt08l4GTKlAk5c+bUFnUREfGooAtenIFL9uzZkSZNGj3RuikgvHTpEsLDw839XLlyuePDioiI3FayYEsVOQOXLFmyePt0Akrq1KnN/xnA8PoqhSQiIp4SVAW7zhoXrriI+zmvq2qJRETEk4IqeHHS2ABdVxER8V9BGbyIiIiI/1LwIiIiIn5FwYuIiIj4FZ8MXhYuXIjChQvjkUcewejRo719OiIiIuJDfG6r9PXr19GzZ0+sWLHCDFYMDQ3Fs88+i8yZM3v71Pyq7wq3hSdLFvPbe/XqVaRIkSLeHy+h7yciIoFl715uegEefti75+FzKy/r1q3DY489hjx58pix2PXr18eSJUs8+kR/8eJFr9z4ueNznp999hkeeugh01OlVKlSmDVrlnnbzz//bHZQ8TqVK1cOKVOmxOrVq1G9enW88sorJhjMmjUrateubR6/cuVKVKhQwTyODeV69+5tgkanO72fiIgEr2nTgNBQoHlzICIiwFZeVq1ahc8//xwbN27EkSNHMHfuXDRu3DjGY4YPH24ew7czUPnyyy9RtWpV87bDhw+bwMUpb968OHToEDyFnWHTpUsHb7hw4QLSpk0bp8e+9957mDNnDkaMGGHSabzO7dq1Q7Zs2VyP6dWrFwYNGmQCHLbqp/Hjx+Oll17CmjVrTADEa8mAsFOnTpgwYQJ27tyJrl27IlWqVOjfv7/rY8V+PxERCU6XLwM9egAjR9r7fMo8fx5IlSqAgheuKHBVoHPnzmjatOktb58+fTp69OhhApgqVargm2++Qb169bB9+3bkz5//tk+Uwd6Xhdd08ODBWL58OSpVqmSOMUD55ZdfzPV74YUXzLEPP/zwllWSQoUKmRUbp759+yJfvnwYOnSoua5FihQxAeM777yDDz74AEmSJLnt+4mISPDZuRNo0QL480+bLurbF+jXD4hVlZDo3P7pGYjwdid8Eu7SpQuef/55c5+rLkx3cEVh4MCBZtUl+krLP//8g4oVK97x4125csXcnM6dOxfvrrBcAfGGuHb6ZWAXERFxS2DCWpQyZcq47jNlFFvsYzt27DABUPSAkEEkrwGvNQPIO30sEREJHhMnAi+9xBfQQPbswOTJQK1a8AmJGjvxyZbpJNZYRBcWFoa1a9eaf7MWY+vWrSaAYcHu4sWLzYrAnTDgGTBgQILPiU/icU3deMvNmzfN/xctWhQjpUasW/nrr7/Mv2/3dcQ+xpWt2CtZztWu6Md9/ZqIiIhnMFh55RVg3Dh7v2ZNYNIkDt2Fz0jU4OXEiRNmF0yOHDliHOd9Tns2J5QsGb744gvUqFHDPGmzjuNuQxT79OljCkujr7wwLRJIihUrZoKUAwcOoFq1are83Rm8xPVjzZ49O0YQw8CRxdGxAyMREQku27bZNNH27QCrCJgiYqooaVL4FK9krW73yj/6sYYNG5pbXPBJnbdAxsDirbfewhtvvGECuieeeMIEaQw6WGxcoECBOH+sl19+2aTqXn31VbOjaNeuXejXr58JAJ31LiIiElwcDmDsWLviwgJdrrJMmcLdp/BJiRq8cNtt0qRJXassTuHh4besxsTXsGHDzI0rO4Hoo48+Qvbs2U2abN++fWY3EXvgvPvuu660UlxwdYWpuLffftsUVrN/DmuQuJtJRESC04QJQJcu9t9hYbbehXUuvirE4cF9sFxNib1VmsW3ZcuWNbuNoqcyGjVqZJ6Y7xdXJDJmzIizZ8+ampnoWPS6f/9+FCxY0GwNFvfS9RUR8U8REQA7ljz7LPDOOzZllNju9vzt8ZUX7lrZyxZ8kRgsbNmyxbzC504Wpifat29vdrNw18vIkSNNLUe3bt3cfSoiIiJyG1y2mD8faNDA1rPw9Tz3zSRPDr/g9uBlw4YNptjWyVlM27FjR4wbNw4tW7bEyZMnTU8SNqkrXry4SWPEp25DREREEo7twTg6kL1JWZRL/hK4eCR4YWv5e2WiWDTKmzsFes2LiIiIu1SrZrdCp0/vn9fUozUv3qCaF+9RzYuIiG9yODh+h5s2oo6xy4a3BywmtOZFe2NFREQC2JkzdphihQrA8eNRx30pcIkvBS8iIiIBav16Owl69mwbuEQ2s/d7ARO8sN6FW67Lly/v7VMRERHxeproyy85u467foGCBYE1a4BGjQLjGxMwwUv37t3NAMP1DDNFRESC1KlTANurvfEGcO0a0LQpsGkTEEiv7QMmeBEREQl2//d/QOnSwIIFQIoUzEoAM2cCmTIhoCh4ERER8XM3bwKffWa75B48CBQqBPz6K1uTsNs9Ao6ClwDFRoCchfT333/H6/2aNWuGwYMHe+y8RETEvY4fB555xrb1Z6uzVq2AjRuBMmUC90oHTPASTAW7169fNzOj7taQj3OiGjRogAcffDDGcU6i5sypunXr3vb9PvjgA3zyySdmv72IiPi2f/6xaaLvv7ct/keOtNOg79Emxe8FTPASLAW7u3fvRvPmzc1Yhc2bN6NDhw4mQClZsiQ2sSILHGd+GWPGjMHzzz9/y/t/++23aN26NVasWGFmSsXGj8OPN3ny5ET5ekREJOHy5AHKlQMKFwZ++w3o2jUw00QBG7wEA86Cql27Nj7++GOcOnUKEydOxPjx47Ft2zakTJnSDLyk77//HsmSJTODL6O7ePEipk+fjh49eqBmzZpm1tTtNGzYEFOnTk2Ur0lEROLn2DHg/Hn7bwYq48dzriBffAbPlVTw4nTx4p1vnBUe18devhy3xyYAg47HH38cjz32GCZMmICjR49iwIABSJs2Ldq2bWtWnljrsmrVKjO1OzYGLjlz5kSFChXM48eOHXvbOVR8+7p163DlypUEnaeIiHjG6tU2TdStm+3lQtxJlC5dcF1xBS9O/M7f6cZN8tFlz37nx9arF/OxrDm53ePi6fDhw5g1a5aZzk2c3M0VFmdNy02WmgNInTq1KdLNnTv3LR+DqSQGLdS4cWOEh4dj2bJltzwuT548JnBhcCQiIr4jSRJboLtli237H6wUvPiJlStXmgCFqy7OXUEsqu3UqZO5/8cff5hi5TRp0pial1Ss3Ipm165dpli3TZs25n66dOnQqFEjUwMTGwMgunTpUiJ8ZSIicjfXr0f9mx1zv/vOtv1/4IHgvW7JEEC7jXi72w6cu7pw4c5vS5o05v3w8LuHxdHFc6vynRw6dMj8n2mf2BjEzJs3z3z9lDVrVpw+ffqWVRcGN48++qjrGFdhmjZtah77QLTfAtbTULZs2dxy7iIikjBLlwIvvAAsXAhEvnZF7AX+YBQwKy/3vdsobdo732KtYtz1sZGrFvd8bDw5g5Z/uC8ultdffx1Vq1Z1raqUKVPGXIvoW6tZI+N8u1OdOnWQPn36W3YWbd26FXnz5jVBkIiIeGe15b33+Hfavgb+6CPf+C44HA7Mnz8fW5i38qKACV4CXb169UzTOe40cta3nD9/Hl27djUrJdOmTTP9W5xBCXcgOVdfFi5ciGPHjpnt1QxMnLedO3eaoIerMtGtXr0aYWFhXvgqRUSEr1Fr1gQ++cQW5b74IjB2rPevy6pVq1ClShVTM/nmm2/edsNHYgmYtFGgy5IlC5YuXWp+YEJDQ/Hwww+b461atTJ9X6IrUaKE2W00Y8YMvPjii67ghNus74Q9YvhxIyIiTAO8JUuWePgrEhGR2Nhsjl0vTp4E0qcHRo0CWrb07nX6/fff0adPH7NJxFkXyZ2vLNNgWw5vUPDiR9hAjgFMXLz//vt46623zMrMd6zuiiMGOhUrVjQ/mCIikjg4/ZlpIs4notBQtrewM4q8Zd++fabr+pQpU8wqCwMVPqfw+SVXrlzeOzEFL4Grfv362LNnjyn0zZcvX5zfL3ny5BgyZIhHz01ERKKw2TnnEXEiNL3yCjBoEJAypXeu0rFjx0yJwjfffINrjKrA1Z+W5lghb0ZT0WjlJYCxkDe+XmBZu4iIJIoFCwB2vGCJYsaMXP2+tbVYYjl37hy++OILc2NHdmL9I2flsazAlyh4ERERSWSsdX3rLWDwYHufM4WZJipYMPG/FVeuXMGIESPMUN4TJ05Enk95fPrpp2aUjC8KmN1GwTRVWkRE/Bs3hzoHKL7xBvDLL4kfuNy4ccO00ShcuDDeeOMNE7iwFxi7uf/2228+G7hQiMObe508tOyVMWNGnD17FhlizQTnTpr9+/ejYMGCt3Sglfun6ysicq+/k1Gtw1hOwllFiR0jOBwO00Lj3XffNW0ziCNl+vfvj86dO3ttB9Hdnr8DduVFRETEV3HO7auv2kAlsgYWyZMnfuDyyy+/mP5eDRs2NIFLpkyZ8O9//xt79+41O4m8FbjEl3+cpYiIiB/jnNuJE4GzZ4Gffkr8Fv9//vmnWWnhiouzVws3dfTq1SvGeBh/oeBFRETEwwoUACZMALiwkZiBy99//216tUyaNMmki5ImTYouXbqgX79+JlUUb+fPM78D5MkDb1LaSERExM0uXwa6dQN++CHqWMOG7MGVOJf6+PHj6NGjhynGnThxoglc2I2do2PYvyVegQtH0nC5iK1/OWfvnXfgbVp5ERERcaNdu4AWLYA//gDmzmWn2gTN402Q8+fPY/DgwRg0aBAuXLhgjj311FOmV0u8d+Pu3g2MH2/zXQcPRh3n4F8GNEm8t/6h4EVERMRNJk2yKy7s8ZY9u72fGIHLlStXMHLkSHz00Udm1YXKli1rerXUqlUr/h+Qqyw8eadMmWwb4I4dgYoVo/Z5e4nSRuJy8OBBVK9e3fTL4RylmTNn6uqIiMTBpUvAc8/Z53wGLjVqAFu2cCCuZy/fzZs3TT1LkSJF8Nprr5nA5ZFHHsH06dOxbt26uAUuN27Y/BZzXU7FitmVFea52D3vyBFgxAiAc++8HLgE1MoLm9TxxqY7kjDcIvfll1+idOnSCA8PN+2gOSMpbWKtd4qI+KFt22yaiNkUPq9/8AGH4wJJk3ruczocDixevNjsIPqD+SnADEtkIe5zzz1n5tTF6cSZFuIKC4OTqVPt6gq9+KKdW+DlAYwBv/LSvXt3bN++HevXr/f2qfiEkydPInv27KbSPK74g8/Ahfi+mTNnxqlTp8z9Zs2amTyqiIhYbPE6dqxt7c/AhbWsrGvt39+zgcvatWtRrVo1PPPMMyZwYWM31rSwV8uLL75498Dl5Em+2rcnXbw48PnnNnDJksXu43bKnNlnA5eACl6CyfXr1zF37ty7rjLxB7lBgwZ48MEHb/mhDwkJQd26de/6OTZs2GCWI50TqbnVjnMv2AFRRCTYsRaW5R9MFTHbwvQQ00SebDrHnUKNGjVClSpVsHr1atMp/u2338a+ffvQu3dvpEmT5u4fIDycrXTt2OoNG+y+7UaNgDlzgMOH7WqLn1Dw4md2795ttrsVL14cmzdvRocOHUyAwhqVTZs2mcdcvnwZY8aMwfPPP3/L+3/77bdo3bo1VqxYgQOcw36HVRt+XBZ/OfHj8/NMnjzZg1+diIjvY5amXDm7CYdlIR9/bEtGcuTwzOfj32q27eff4QULFiBJkiTm7/uePXvw2WefmVXy2/r9d/7Rj7rPCmKuuHCF/csvgUOHgHnzgCZNgBQp4FccAebs2bOc1WT+H9vly5cd27dvN//3R4cPH3bkz5/fsXXrVsevv/7qeO211xw3b950XLhwwVGuXDlHsWLFzONmz57tyJo16y3vz8elS5fO8dtvvznq1q3rGDBgwC2PiYiIcFStWtUxYcKEW97Wv39/87Y78ffrKyISFyNHMmHkcOTO7XCsXOm5a3b8+HHHG2+84UiRIoV5XuPt2WefdezYsePO73TsmMPxn/84HKVK2ZNMntzhOHEi6u23eW70h+fv2LTyEonV4fG9Xb8eFQTy3zwWvVj7bh83Idhw6PHHH8djjz1mJoEePXoUAwYMMAW1bdu2NTU/XDVZtWoVyvFlQSysPs+ZMycqVKhgHj927FhT9BUtkEWnTp3MJNH2LJmPhe/H6nVuyRMRCVZc1B40yKaJnnzS/R+f/Vk+/vhjPPzww/jPf/6Dq1evokaNGmbS8+zZs83OohiuXrWpH6aA2PmWY6q56sLVFB6LXstyj4GH/kLBS6R06eJ/Y/MhJ/6bx2K3fWbJye3eN74OHz5sxpR3ZJIV3IZXA99//72rpoX1Kc55FSzSvV33RKaSGLRQ48aNzY6iZcuWud6+Zs0aE+DMmzfPFO7yxnkYTnny5DGBC4MmEZFgsXmzrWk5fdre546iN98EsmVz7+dhkMJds4UKFcL7779vagzLlCmDH374wfyt5gvI2xo1CmjaFFiwwL6S5uNYlMtCXLa8eOghBJqA2Sod6FauXGkCFK66OHf/8ObEinN2T2TBFmteWMgV3a5du0yxLldbKF26dKbwizUwzj4ATzzxhCsIuh0GRnSJDQ1ERIIA/yS2a2d3E737rm114v7PcRPTpk0zAQuLb4mrLlx9adGihalxcWFAMmkS8OijdlWFuE+bu4acTeSKFkWgU/ASKbKLcrykTBn1b9Y78WPE7pYcj53Kd3WIhVXgVryct7yN0TlXSxixU9asWXHa+RIh2qoLg5tH+QMfiaswTZs2NY+Ny1RR57bpbO5+uSEi4qP4N50DFT/9FPjkE/d+bKbquYLet29fbGEOCiz6zWF2d7IgN4WziDYiApg/Hxg3DvjxRxtRVa0aFbzwb/L+/T7RPC6xKHiJdL992LjjjLfY3NXfzRm0/PPPPyYij45jzatWrYo2bdqY+1xmZMfF6FurWSPDrXTR1alTB+nTpzc7iF7h1rl72Lp1K/LmzWuCIxGRQMV2YZxPxBUXKlvWZl/cmR6aOnUqvvjiC1dqPkOGDOjVq5epbXQ1Bv31VxuwTJsWs26lcmWgQwfbaMYZsARR4EIKXvxEvXr1TOM4LiNyFYXLiBzA1bNnT7MiwiVH9m9xBiV9+vRxragsXLgQx44dM9urGYBEx6CHHy8uwQv7CoSFhXnsaxQR8SbGAl99BfTqZWMBZunLlHHfx+ffZE50HjJkiKljJAYq3bp1M3+zs7BRXHRs1bt0qf03e24xJcSg5ZFHEOwUvPgJ/lAvXboUb775pmnb71x9adWqlen7El2JEiXMbqMZM2aYbosMTqj2XYZssEcMP+6dREREmMZ4S5YscdvXJCLiK5gV79zZ1rwS618LFnTPx+YmCo5eGT16NC5GbjflpgrOInrhhRfwANND3PXB/BTrErljiNg0jqvubNNfvbpXpzj7HEeACeQ+L/GxaNEiR9GiRR03btxwy8cbOnSoo3bt2nd9TDBdXxEJHGvXOhz589u2KClS8O+dw3Hz5v1/3HXr1jlatGjhSJIkiatPS4kSJRzjx493XOHfSTaJ6dzZ4UiXzn5y3j791BGszsajz0vArLxoMGNMHKjI7oss9HW2+L8fnJXBpU4RkUDButcvvrC7iLjDmAvaM2YAd1mEjsPHvIlFixZh0KBBpueWE1e+33rrLdQuVw4h/FvKqc0ssnXidmamhFq2vM+vKjiEMIJBAOHOGw6pOnv2rCmAip362L9/PwoWLHjLVmK5f7q+IuIvTpywJSSLF9v7jBk4ESWhPdzYomLixIlmgC1bUzhf9HEcS8833kCpyKG3pvCWqSDuIEqfHmDanyfC3UNBVnQbn+fv2AJm5UVERCQuVq8GWre2o33Y8uK//wW6dk1Y7HDixAkMHz4cQ4cOxfHjx80xPgG/9OKLeDM0FFkXLeKWUDbrsu+QMaMdO503r+2xca9hinJbCl5ERCRo0kTs18JNPDduAIUL2zRRyZIJG5LL1v3jxo0zq86UP39+DGjXDq2vXkVKbm/+7LOod2CXO6aK6J133PUlBS0FLyIiEhQ2bQLee89WxrKHC7vlxmdcC6ssOEaF/Vnmz5/vmg1XtmxZDKpTB9VWrEDIv/4V9Q5cZWHXW+4WCoKut4lJwYuIiAQFzqv9+GNbcsJt0XFNE924ccO0imARLocjEjctP1u3Ll7t3RtPPvkkQrjHmoELtzPXqWPrWNgBV/WVHhGUwUuA1Sj7DF1XEfElTA0xTcRi3EKF7DHuLIrPdGfOg2N6iJs9qFTy5Pi0aFE8dfQokpcvD1SrZh9cvz4weLD9ZLcZjCvuFVQdb1j5TRos6BnO6+q8ziIi3tSnj00TcW7htWtxf78jR46YeUOsYWEjuTP79+PttGlxMHdubLl2DXX/+APJw8OjtioR/+698YYCl0QSVCsvSZMmRaZMmRDOHzqwyDuNq6W+3N+KCwMXXldeX15nERFve+01O5OoZ08bW9zLtm3bTD0L571x/hBNz5ABTS9dQlJ2xuWNQ+yeftqmhfh/8YqgCl6iDzh0BjDiPgxcbjf1WkQkMbDRHCeYOGMK7kZmyxXncOY7vfhavny5CVo44ZlltQxbqlSpYsaxNFq8GElGj7ZDjhiwcAAupziLVwVVk7rYBVjX4rOOKHfFVJFWXETEW9izhb1b2MOFY4IaN7774/n3n/PfWIT7z5YtaAOgEwDOYfxz1CiUeP55+8C9e5kTT9h+aokXNamLAz7R6slWRMT/ff+97azPrrlsWstC3TvhC1sOSBz2n/+g1KFD6M9aW74Acz4gRQqUiP6a3lnpKz4l6NJGIiISGLh4zoJcZy84ZnamTwceeeTWxx48eBBfffUVRo4ciZznz2MdgKzRH1Chgk0LsS9L5syJ9SVIAil4ERERv3PggE0TrV1r73fvDgwadGtblc2bN2PMxx9j57x5WMYWu+yEW6QIUh07BkeqVAjhkg2DFjWR8ysKXkRExK98951tWnvqlG1iO2YM0LRp1NtZyvnjggVY17cvym3bhq8AcOpQu2rV0LNXL9StWxdJ9u0DChZkDYE3vxRJIAUvIiLiF7h7mb1b2AuO2COOI4QeesjevxIRgR8//BCXvv4adU6fRp1o75uuVCn8NGmS3YJEqmXxawpeRETE57HBLctR1rFYBUCPHsC//223QZ86dQpff/01kg0ciF4XLrje53T69EjSqRMyvvoq0t2uEEb8lk8GL02aNMHPP/+Mp556CrNmzfL26YiIiBdduQJUrWq3Q2fKBIwbBzSqfQnhQ7/G2N9+w4cLF5pGmaUAvBISgv1lyqBAv3544Jln7KwhCTg++V1lO+YJEyZ4+zRERMQHpEwJDBwIPP64A1u++Q2VxjTEpYwZkf3NN5FzxgwTuJQuXRpvTZyIZCdO4LGNG5GuYUMFLgHMJ1deatSoYVZeREQkOLE33MmTQMWKNmfUatcYNNg7GplaHnM9Zh+3SxcujJ+GDUPNmjU17iWIxHvlZdWqVWjQoAFy585tflDmzZt3y2OGDx+OggULIlWqVChbtixWs+WhiIhIHPz6KxAaCjz7LPDbb/sQXrYskn/yCTKdOIbzAMaFhODTevVw6fff8fzOnabEQHPqgku8V14uXryIUqVKoXPnzmgafW9apOnTp6NHjx4mgOFsiG+++Qb16tXD9u3bzYROYkBzhUnMWH788UcTFImISBBiH5aff8Zj42YiU4Z/4/zFA3j88TD0xmnUBLAoSxbkffVVtO/WDTly5PD22Yo/BS8MRHi7k8GDB6NLly54PnIuxJdffoklS5ZgxIgRGMikJYCNGzfCXRgERQ+EOBtBRET8yJ49+N9/5yP3nK+Q/PA/SA+gFBZgIY4hJOQmfq9TF1W6d8cX9epprIu4v2CXI8QZmISFhcU4zvtrnW0Q3YwBEQcxOm/58uXzyOcRERE3On8eGDUKjipVMOnRAXhsaDd8ebgVzgD4GsCZTJfQq9eb2Lt3LxZ//z2eeeYZBS7imYLdEydOmGnNsZfzeP/o0aNx/jh16tTBpk2bTIoqb968mDt3LsqzG9Ft9OnTBz179oyx8qIARkTEt53fvRtJX3gdr2IIvkUXc+wL1Maiimvw/CsvY2mzZqZuUiTRdhvFLpxiq+b4FFMxzRRXKVOmNDcREfFR27cD48ezaBKbu3QxZQQTJ25EFqzDIRRnsQvKll2MkSPzITTUM6v0EljcGrxkzZrVLOvFXmUJDw/3eHHVsGHDzI0rPyIi4mUcPDR1qg1a1q83h66GhOCpYcNwGh0BrMYhpEGGDBfBrv0NGjzj7TOWYK15SZEihdlJtHTp0hjHeb9y5crwpO7du5sdTesjf0lERMQLVq0CmjUDcuUCXnnFBC7XAMwH0MyRBudC2IB0HIA0qFXLgd2706JBg7T6VolnV14uXLhgCqic9u/fjy1btiBz5sxmKzTrT9q3b49y5cqhUqVKGDlyJA4cOIBu3brF91OJiIi/bHGObMN/c8MGJJk92/x7M4DxAKYASJqzNm7cmIIbx7Oah374IWsWQ9S9XxIneNmwYYPpgOvkLJbt2LEjxo0bh5YtW+LkyZP48MMPceTIERQvXhyLFy9GgQIFEnaGIiLie44fB6ZMsYOGXnsNR+vVw5gxYzBr+HC0jwxa/gwJQVhYHbQo9CnGjCmJiIgQsJUXs0lPPuntL0D8WYiD1bQBIHrNy+7du3H27FlkyJDB26clIhI4rl4FFi2yAcvixcD16+bwnzlyIPTkSVyPvJ8lSxY899xzePHFFzF8+MMYPNi+e926AMfWZcvmzS9CfBV3C7PlSVyevwMmeEnIFy8iInFMC3GVnZW1HDgUaWuqVBgREYFprM8FTKnASy+9hObNm7u2Oa9ZA9SqBQwYALz1lmYlinuev31yMKOIiHjZmTNApkz23yxS2brVBC5n0qTBt9euYfS1a9gREYE0adKgXbt2JmjhZGe+HN61CyhSxL5rlSpmriJy5vTqVyMBRsGLiIhYERHAggU2LbRsmYk6IjJnxqxZs7Dm6FHsB/DTpUtgQ4qiRYtiyMsvmw0afLVMly6x/tFmlrjx87HH7IdV4CLupuBFRCSYcankt99sP5Zp0+yKS6TpXbui+2+/mU0YlCxZMjR99lm8/PLLePLJJ29pPspM0dmzthRm06ao4EXE3QImeFGTOhGReNqyBWjVyuZ5Il3OmhXzM2ZEv7/+wm4W5QJm5MoLL7xgBu7mjLWMwtjn2jX2+bLZJZbF/P03UKGCvhviOSrYFREJFszrHDoEPPKIvc9Vlpw54UiSBFsLF8bAw4cxLTwcjmhz5ljL8vTTT5tVl9s10X3uObt7aNSoxP1SJPCoYFdERKKWRn75xaaFZswAChc2BSncaPrLn3/i58qV8Z/Vq3GaqzCAaTjq3OZcqFChO17FX38FWrYEDhywqy7vvAPc5eEibhUwaSMREYmGuRs2VWHQsm+f6/DN48cxdtAg/GfcOGzbts11vGLFiqaWhducU6dOfddd0+zb0qePrW15+GFg+nQFLpK4FLyIiASa994DPvkk6n66dDj11FMY63BgwLJlOP/22+Ywtzm3adPGpIZCQ0Pv+WFPnAA6dbK7iYgrLyNHAmqpJYktYIIXFeyKSFDiUsjPPwOPPgrkzWuPlS0LhITgZvXq+LVwYby/eTOWz+doRKtIkSImYOnQoQMyOXu53AMzT61bA//8A6RMCXz1FfDCC+bTiCQ6FeyKiPgjDshlSoipIRae9OsH9O9v3rR/1y5M++orDJ45Eye4XBK5zblJkyYmNVStWrVbtjnfLTb697+B998HbtywMRJLZ0qV8uhXJ0HonDrsiogEIDZRYeTAoIV9950yZsT1mzcxb9YsjB49Gj/++KMpyKU8efKY4ltuc86VK1e8Pl14ONC+PfDjj/Z+u3bAiBEmCyXiVQGTNhIRCWisjuVOoWPH7H02ValTB4dq1cLQgwcx5uuvcZyTniPVrl3brLI888wzt93mfC/sW9ekCXDkCMD63aFDgc6dlSYS36DgRUTEF+3YAXz3HcDiWqZ4GIA884zZo3y1dWvMT5cOQ2bPxuo333S9C1dWOnfubLY6P8xtQPche3bbFqZYMbvYo2654ksUvIiI+Ap2fWOLfqaF1q2zx2rUAMqXN//c0qULRqVIgcmff24m71KSJElMEzmmherXr5+gVRYnBitp0th/Fyxo00UMWtKmdcPXJuJGARO8aLeRiPhtOmjJEjsMkUMRr161xxmE1K+PCxERmPT116aWZePGja53K1iwILp06YJOnTqZupb79dNPtr6Fp1Gnjj2mFv/iq7TbSETEm5YuBcLCou6XLg1Hhw5Y/8gjGDF7NmbMmIFLXBIBO9mmMDuGuMpSs2ZNs+riLt27A8OHAzVr2kBGW6AlsWm3kYiIL2JB7ZQpQNKkwCuv2GOMFsqUAapXx+mGDTF20yaMHjUKO1jzEqlo0aLo2rUr2rdvj6xZs3rk1L74wow5wltvKXAR36eVFxERT2IaiC1pmY/hlGamiRglHDxoUkM3b97Esp9+wugxYzB37lxc44jmyO63LVu2NKsslSpVinNflrj64Qdg7NioWErE27TyIiLibRx0+O23Njo4eTLqeLlypsf+oQMHMHbKFIwZMwZ/cw5RpLJly5pVllatWiFjxoxuPy3GRmw4x8ZzVK0a8PLLbv80Ih4VMAW7IiI+hYHLkCH232wO1749rrdti8V//41Ro0Zh8WuvmVUXYpDSrl07U4BbhikkD+FiT6tWwNq19j6Dluee89inE/EYBS8iIvcjIsLuEmJaqHdv4Mkn7XF2dGOL2k6d8FfBghgzfjzG1a2LI+z6FunJJ580aaGmTZuaNJEnsWUMhypyNzYHKY4ZAzRr5tFPKeIxARO8aKu0iCQatt5nC1oGLNOnA2fO2OM5criCl4iiRTGvcWOM+vxzLF++3PWu2bJlM9ubucpSmB1zE6Hkpk8fYPDgqKwVT/mhhzz+qUU8RgW7IiJxdeWKjQIYtOzeHXU8Xz7bJKVjR2y7ds30ZJkwYQJOcZmDf2hDQhAWFmZqWRo0aGC2PCcGltK0bBnV765HD+DTT+1UaBFfo4JdERF3YV2Ks58Kgw7mW/76yw78adrU5GIulC+PGbNmYVTHjvj1119d75ovXz7Tqp8t+wsUKJCo35O5c209CxeFMmWy8VajRol6CiIeEzBpIxERt6aFfvnFtulnx7Zdu+xyBbcr9+9vtuw4mjbFhl27zCrLlMaNceHCBfOubM/P1RWusnC1JWki70Pm4hDHITlrhR9/3E4cSOTYScSjFLyIiETPs0yYYIOWfftiNkWJXLY4/fTTmDx5MkY/+SR+//1310MeeeQRU3zboUMH5GQfFy/WD7OtDDGI+eQTIHlyr52OiEcoeBER4cwgtpb9+eeoa5EuHdC8uUkLOapUwaqVK80qy6xZsxDBCAFcjEmJZs2amVUW7hxydyO5hGBrGE6BPnoUePppb5+NiGcoeBGR4Kxj4VTmBx6w9zk2mYELgw+26+/YEXj2WRy7cAHjx4/H6Oefx549e1zvXrJkSROwtG3bFg84P4aXMI7q2dNOf+Z8Iipb1qunJOJxCl5EJHjs3WtTQkwNVaxolyioSBFg1CgzIPFGnjz48ccfMbpDByxYsADX2c7fLMSkQ+vWrU3QUq5cOZ9YZSFuex4xAkiVyvZt4W5tkUCn4EVEAhtXWBikMGhZsyZmZStvkfuGD4SF4dtvvzW3g2xFG+nxxx83AUuLFi1MAONrOnQAVq+2GS4FLhIsAiZ4UZM6EbkFh/gMGmRzK8Qtz2FhNi3UqBGuhoTgu9mzTS3LkiVL4OAuIwCZM2c2E5xZgFu8eHGfurCXLgH/+hfwzjtA+vQ20zV6tLfPSiRxBUzw0r17d3NzNrkRkSC0YweQP7+tYSHWozBwKVbMBizt2gG5c2PXrl0Y06+fqWcJZwv/SDVr1jSrLI0bN0Yq5mF88Mtr0QLYuhX43/+AiRO9fUYi3hEwwYuIBCl2sWUjE6aF2EqW3dgYqDhzKlWrmp74lyMizE4hrrKsWrXK9e7c1swmcmzX//DDD8NX8cvjIEWuvHAnNkcniQQrBS8i4n9YRLtkiQ1UOBSRA3yIDeHY/TaSI0sWrN+3D9++9BKmTZuGs6x/MdmjJKhfv75ZZeH/2VjOV128aIMW1hhTrVrApEmqb5Hg5ru/sSIit8NOto8+CkSbzoxSpexqS5s25lmdqaBJkyaZ4ttt27a5Hvbggw+aFRautOTJk8fnr++ff9o00c6dtlznww/t4OpEbtor4nMUvIiIbzt+3O4SatzY3ueOH25t5upL27Y2aCld2mxp/uGHH0zA8t1337m2OLN2hY3kOGOoWrVqZtXF17FumCOUXn3Vluzkzg1MneoaWC0S9BS8iIjvYRqIPe5Z6MH/89n8n39ssQexUjV7dtP3nsW3Y3v3NsW3R9lWNlKFChVMwNKqVSu/KuI/fx7o1g2YMsXer1vXpoyyZfP2mYn4DgUvIuIbGKBs3mzrWPjMffJk1NvKlbNposjg5XyGDJgxYYJZZVm7dq3rYdmyZTNbnJkW8rUtznGxZYtNE7GZL1ND3BLNqQV+sFgkkqgUvIiIb2BehGkgp1y57NZmpoUee8z0YPll9WoTsMyYMQOXuO3G1OgmNUW3DFiefvpppEiRAv6qVy8buOTLZzdQVa7s7TMS8U0KXkQk8bGQ47vvbD+W+vXtsXr17FRB5kkYsNSuDSRLhkOHDmHCwIEmaNnL9v6RChcubNJCXGnJxUAnAHz7LdC3LzB4MJAli7fPRsR3hTicLSUDhLNJHbdEZsiQwdunIyJO/FPj7MPCZYUzZ1iYAvz2W8ygJlUqXLlyxRTdMmBh59ubHKQYOV+oZcuWJmipVKmSz8wXup9h1kuX2h1EIsHuXDyev7XyIiKexUJbNiZh0LJrV9TxvHmBp54Cbtxw7f39Y/duE7Bwm/PJaDUvVatWNQELdw354nyhhOD4pCpV7Hglluc884y3z0jEfyh4ERHPYoc1pogodWrg2WeBTp2AGjVM0HL69GlMnTrVBC0buRQRKXfu3OjUqZO5PfLIIwH3XWJdyyuvAPv2AU884e2zEfEvCl5ExH1pIfZj4fbmd98FCha0x1m/whQRA5ZmzYAMGUwaaNmyZRg7dizmzJlj0kSUPHlyNGrUyKyy1K5d26c73yYEM2Qsz+H4Jfr0U7vo5OfZL5FEFzB/GTRVWsRLOCGQjUgYtDhb83NZ4YMP7L+bNrU3APv378e4L77AuHHjcODAAdeHKFGihOl827ZtW2TNmhWBGNexCJe1LeXLAytXmhY1rEcWkQQImF8dTZUWSUQsrJ0+3QYsK1ZEHWc9SvPmQJ06rkOXL182qytMCy1fvtx1PFOmTGjTpo1ZZQkNDfX74ts7YekOF50WLoyK67jQxOBFRII8eBGRRMQiWxZscM4Qg46aNW16iPUsadOaniwb1q83AQvrWZwDERmg1KpVy/Rkady4MVKzBiaAMYvWqpWtWU6ZEvjyS+DFF5UmErlfCl5E5O7YW4VpIRbTcvmAwQr7s/TsCbAhXPv2riIOMxDxm29uOxCRAUvHjh1RoECBgL/i3Nn92WfAe+/ZOI9zJGfMsPMjReT+KXgRkVtxpWTmTJsW+uWXqOObNgFly9p/Dxhg/mcGIi5ceNuBiE2bNjVpoerVq/vFQER3zZHs0AH44Qd7n02DR4wA0qf39pmJBA4FLyISc7jO558Dc+eyWMUeY9ARFmbTQsWKuR5qBiKOHXvLQMTy5cu7BiKyriWYrFoFtG4NHD5sd4UPGQI895zSRCLupuBFJNhFaxJnhh86xxkzUGHAwvlCuXObQ+fPn8fMb781qyxrWNARiTuEOnTo4LcDEd1xCQcOBPr1symjokVtmigIL4VIolDwIhKMTp2yu4XY9bZ6deDf/7bHOU/o7bftjiFOcg4JiTEQcebMmbh48aJ5KNNAHIjIVRZ/H4h4v/7zH+D99+2/ubNo6FBbFiQinqHZRiLBgrUoS5bYgGXBAuDqVXucxbZ//31LbsMMRJww4ZaBiI8++qhrICK74ArAeI4x4Kuv2noXEYk/zTYSkZg+/tguBxw7FnWsZEm7TNCmjStwcQ5EZC3LDz/84BqImDZtWtdAxMqVKwdsT5b4pIkmT7YZNZYEcZWF3XODpCZZxOuUNhIJRCdOAFmyRK2msKCWgUu2bHb7C2tZSpd2PfyPP/647UDEJ554wgQszZs3D5iBiO7olsshitxNxEvKLBspcBFJPApeRAIF00CLFtntzfw/u9lWrWrfxoZy3DFUr56rteudBiLmypXLNRCRKSKJifFgkyZ2Bzm75YpI4lPwIuLvywCbN9s6Fu4SirZqEiN4KVLE3G7cuIGlP/xg0kLz5s3D1ci6Fw5EbNiwoVllCQsLC7iBiPfr2jXbJdc5a7JrV7v6opIfEe/QXygRf8VtzVxN2bo16hhHFrMQg2mhxx6L0ZOFwxBZgHuYTUiiDUTk9uZ27dohG1NKcouDB22L/0OHbJz4wAN29UWBi4j3KHgR8adhiDt3RtWq5Mhht7lwaE7jxjZg4VbnyFUTzhOaMWOGWWX5v//7P9eHyZw5s5nezLRQmTJlgr749m44DYGXlTvLM2SwcaJzMUtEvEfBi4ivp4XWrbNpoWnTbFUoV04YsPDfbOH/0EN2OcDM1LmJFcuWmYCFk5w50dnZk6VevXomYGnQoAFS8v3ljphN69MHGDzY3udEBLbFefhhXTQRX6DgRcQXscBi0iQbtOzaFXU8b17gr7+i2vRHzhnat2+fSQuxVf+BAwdcDy9atKgrLcRCXLk3trxhmohbn+n1120PP8V7Ir7D54KXgwcPmuZXnE7LosH333/fbNMUCRqjRgEvvmhXXYhDcp591vZkqVHD1cr/woULmDVrllllWcWhOpEyZsyI1q1bm1WWChUqKC0UD/PmAZ07A2fOABzLNHaszciJiG/xueCFAcuXX36J0qVLmwAmNDTUtCBnkyyRgMMAhTOCOHK4VCl7rEoVe5zFFQxYmjWzBRfm4Q6sXrXKBCzRW/WzaVzt2rVNwNK4cWOkZsAjcXblCtCrF/Df/9r7FSvaLN2DD+oiivginwteuLTtXN7Onj27KS48deqUghcJLP/7HzBhgr2x9X7LlvbZkpgSYuonWhMRpoKYEmJqiCkip0KFCpmAhUMR86npSIIwC8fL72x189ZbwL/+5WqHIyI+KN7NrLk8zYI/zjThqz32ioht+PDhKFiwIFKlSoWyZcti9erVCTq5DRs2mAJE/VGWgHDhgm0gV7OmfUn/wQc2cOGqYubMUWkiypcPly5dwuTJk82KyoMPPogPPvjABC7sdMt+LPy92r17N/r27avfkQTavh0IDbWBC78F3F30+ecKXEQCbuWFy9SlSpUyRYBNmza95e3Tp09Hjx49TABTpUoVfPPNN2aXw/bt25GfA+BMjWFZM0Mlth9//NE16I0tyvlqcvTo0Xc9H36c6B+Lg51EfFKtWlFVoGwUwvoVpoVYzxKZFmVa6NdffzUrLNOmTYvx81y9enXX753SqO7B3n1MEV26BEydqo65In7DcR/47nPnzo1xrEKFCo5u3brFOFakSBFH79694/xxIyIiHFWrVnVMmDDhno/t16+fOY/Yt7Nnz8bjKxFxsz17+MPpcFy4EHXs888djkKFHI6PPnI4/v47xsMPHTrkGDhwoKNw4cIxfo4ffPBB8zO+b98+fYvc+K25eDHq/unTDse1a7q8It7G5+24Pn+7teaFrcY5I6V3794xjrPd+Nq1a+MaTJkcfs2aNc2uo3vp06cPevbs6brPV6pKM4lXnD1r+64wNcTBN8TGIM6f49deA9580zUsMSIiAgsWLDCrLEuWLHFNcE6TJg2aNWtmfg+qVatmerSIe8yZY5vOcSs0N3URdxWJiH9xa/By4sQJMzslBzt/RsP7RznVNg7WrFljUk8lS5Z01dNMnDjRtDG/HTbbUsMt8ZobN4Bly2zAwmdGdsElBhzsdsu+LE4pUpjgfOOGDSZgmTJlihmOGH2CMwMWtgbIELm7SNwrY0bblHjPHvutSpVKV1jEH3lktxELeaPjH+zYx+6Ef8Cdr0BF/GLwTZ06UfeLFrUv7TlfKE8e1+Fjx46Z4ltucd4abRZR3rx5TW0Xg5ZHHnkksc8+KEQPUp56CvjpJ+DJJ11TFETED7n11zdr1qxImjTpLass7NcSezXG3YYNG2ZuXPkR8QiuknA7M7c5f/qpPcZdQyy4zZnTFt+WK+dKCzGNunjxYhOw8P/Xr183x7lS2KRJE1N8+9RTT5nfGfEMLoixzT97+BUqZI9xs5eI+De3Bi8pUqQwO4mWLl1q/jg78X6jRo3gSd27dzc31ryww6iIWzDgWLLEPgvOn2+H3jDYeOMNOxiRZs+O8S6///67SQtNmjTJpFKdKlasaFZYWrVqhUwqtPAopoa6d7ffNhoyBPjqK89+ThHx4eCFLcn3sjdFpP3792PLli2mmRy3QrN4loW25cqVQ6VKlTBy5EjTYKtbt27uPncRz+E8IVZ0Tp4MRF9JZO0VV1hiFUtwaz9rWLjKsnnzZtfxnDlzmt8HBi3FnPOIxKOYlWvRAtixw5Ye9e8PvPuuLrpIQInvVqYVK1bcdmtyx44dXY8ZNmyYo0CBAo4UKVI4QkNDHStXrnT44lYrkTsaNoy9AOwta1aH4/XXHY5NmxyOmzddD7l27Zpj4cKFjqZNmzqSJ0/u+l3gv3mMb+NjJHHwWzN6tMORKpX9tuXO7XD8/LOuvoi/iM/zdwj/gwAQveaFXUfPnj2rHRtyb0wDff+9nd7coAHw3HP2+KlTwAsv2G3O9eqZnUJOO3bsMCss3AUXvb6rTJkypo6lTZs2yJIli65+Ijp/HnjpJbtQRqyhnjgRyJZN3wYRf+Es+4jL83fABC8J+eIlSPFHfssWG7BMmcI9/vZ4pUrAHfoRnTlzxnS8ZdCybt26GEXq7dq1M2khdp6WxPf77zZNtHu3LUf6+GM7ZFHtcUQC9/lbmwUluLBqc8wY4M8/o45xpxC3NnOLczRcxVu2bJkJWObOnesaQ8HdQU8//bRZZeHEcxaqi3di0G++AXr0sFOh2VKHm8E4lFtEApuCFwn83ULRG3r8+KMNXFKmBLgDjgFLWFiMx+zZs8fsFpowYQL++ecf1/HixYubgKVt27Ye3/ov9/62tm0LzJhh7z/zjF1IU7ZOJDgETPCiPi8S4yX5+vX22YzPbhwZXKCAfRvb87O2pWVL4IEHYixXzpw506yysMuz0wMPPGBqWBi0hIaGxrnZongWY00GKvz/v/9td67rWyMSPFTzIoHj0CFbpcnmHjt3Rh3//HPgrbdueTg7Oa9YscKsssyePRuXL182xzlLqE6dOiZgadCgAVKph7zPxKTs35IuXVTn3G3bOKXe22cmIu6gmhcJLvv3A+wjxL7vztESqVPbzrfsyVKjRoyHs0/R+PHjze0g2/tHKlKkiCstlCdaa3/xjebGXbqwcJpNL21hLmNKBS4iwSlg0kYSZC/Bw8OjOtxyPyxTPQxcqla1dSzNmwPRqtWdaSGusvzinPhsJgpnQuvWrc1uofLlyyst5KP47Wa5Ene2b9jAbsXePiMR8SYFL+I/OFPImRbiygr3yLLQgXkEHi9ZEnj44TinhRiwNGzYUGkhP1C4sP0W58tnx0eJSHALmOBFBbsBikUOnB3EgGX58qjjadMC3AnEZzOKNkvrXmkh9mXJnTt3on4ZEj8nTwLPP2+3QVerdsu3WESCnAp2xXd9/TXw9tscqGXvc5WF9SusY2E9CwOYSEoLBQ5mAFu3Bhh3PvSQHTMVfbe7iAQmFeyKf+LATwYkuXLZ+1xVYeBSqJCtY2GrfueWZ6WFAg5Llj77DHjvPTYIBB55xO50V+AiIrHp9Yx419mzwMyZNi3EQtrevYGBA6MG1PBlONv2R2viobRQ4Dl+HOjQAfjhB3ufKy/snps+vbfPTER8kYIXSXx8Wb1smQ1Y5syxDTuIw2iOHYv205kMqFzZ/FNpocC1ciXQpg1w+LDd/jxkiN0WraZzInInCl4k8bc5h4YCf/wRdaxoUZsW4nyhaP1VtFso8GPYf/0L6N/fpoyKFLFpohIlvH1mIuLrAiZ40W4jH+4uNn++DU74Upo3Ts5jNSZzAyy+5d5XpYWCytGjNlblAhzxx2PYsBg12CIid6TdRuKZqXlLlti0EAMXdhZjPYtz3C/3wbI3C4cjRlJaKHgwYOFQRWYI06QBhg+/ZaC3iAQh7TYS7+C0ZgYskyfbl9ZOzANcuhR1P3L0r9JCwftjwsCleHFg+nSgWDFvn5GI+JuASRuJl61bF7Nne9as9uU1X1KXLn3btNCECRNw4MAB13E1kQvsUifnj8DrrwMpUtiMIVdeRETiS8GLxB/TQN9/b9M/zz1nj7FuhY05+HKaAUu9evYZKpLSQsGLGcSPPwYWL7ZbnxnEvPyyt89KRPyZgheJ+0vnLVuAceOAKVOAEyds+odVlwxSuM152zYgeXLXuygtJFeuAC+8AHCB7dNPgU8+0TURkfun4EXujsUJrGFh0MJiBSdOdGbHWw47dK6wRAYuSguJE2uyGevy9v77ui4i4h4BE7xoq7SHDBpkb8QgpXFjmxYKC4vRt11pIXFatMg2TmbjOeImM+dGMxERd9BWaYlKC61fb3cLNWtmByDS1q12vC+rK1u2BB54wHXFlBaS6K5dA95918a6qVMDmzbZxnMiInGhrdISd4cOAZMm2bTQzp1RjeWcwQsLcH/9Nca7KC0ksf3vfza2/e03e5/xbsGCuk4i4hkBkzaSeGAvdjbYYMDy00/2PvHlcpMmdrBMLEoLyZ3Mmwd07gycOQNkzAh8+y3w7LO6XiLiOQpeghH3qn70EbBjh73/xBM2LdS8OZAhg+thSgvJvXYSvfMO8NVX9n6FCsC0aVpxERHPU/ASDOv5EycCc+cCq1fbrmAMXt58076tQwegUKEY76K0kNzLvn1AixbAxo32Pn+cOGQxWmsfERGPUfASiC5eBGbPtsW3y5dHHWcAw663FCs1dPbsWcycORPjxo3DmjVrXMczZcqE1q1bo1OnTihfvjxConXKleA0a5b98Tl3Dsic2WYfGzTw9lmJSDBR8BJI/vrLtjKdOdMGME4svmVaqFGjGA+/ceMGli9fbgKWOXPmICIiwhxPkiQJ6tSpYwKWhg0bIlWqVIn9lYgP4o8HV1g4SJEqVwamTgXy5/f2mYlIsFHwEgj7U51dbbkqwpfB9PDDth8LG8k9+GCMd9m1a5drttAh7jaKVKxYMROwtG3bFrlz507UL0N839q1UYFL797Ahx/GaKgsIpJoFLz4I67Xc3WFaSG26Gc6iB56yPZgZwEuXxZHS/GcOXMG06dPN6ssv0bb+vzAAw+gTZs2JmgpW7as0kJyRzVr2jpvjrGqW1cXSkS8J2Ca1EXvsLt7925Tw5Eh2s4Zv3fjhq1fYcAyZ45ty+/svx4eHmOXUNS73MDSpUvNKsvcuXNxhdtDACRNmhT16tUzAcszzzyDlPwYIrFcumSbzvXsqdSQiPhWk7qACV4S8sX7jZEj7Uvef/6JOsbWpaxj4WDEPHliPHzHjh0mYJk4cSIOHz7sOl68eHFXWihnzpyJ+RWIn+FfBf54TZhgF/JWrYqxkCci4nbqsOvv2OGWe07TprX3r161gQtb87dqZZ9VypeP8Wxy+vRpTJs2zaSF1q1b5zqeJUsWV1qoTJkySgtJnAIX/mj172875vbrp8BFRHyLal58xfXrwI8/2rTQ/PnAkCFA1672ba1bA1wp4X7UaCme69ev48cffzQBy/z583GVQQ6/qcmSoX79+iZgefrpp5FCzTckDrhB7Ycforrjsr3/tm1MM+ryiYhvUfDibRx8yICF84WOHo06zl4rzuCFRbkcluh6l60mLTRp0iQcjfY+pUqVMgELV1qyZ8+eqF+G+DcGKWw6t327DWDq1LHHFbiIiC9S8OLNLc4sJoiW4kHWrECbNjYtVLp0jIefPHkSU6dONUHLhg0bor1LVrRr1w4dO3ZE6VjvIxKXFBFnEb36qq0B5w55NmEWEfFlCl4SM1hhoFKlir3PBhlsT5osGfDMM7YnS/36MfqrX7t2DT/88IMJWBYsWGDum29asmRo0KCBCVi4a0hpIUmI8+eBl14CJk+297nawkkS2bLpeoqIb1Pw4mlbttjGcVOmACdO2HlC+fLZt3GiHYtwYz1b/PHHH660UDi3QUdiwS3TQmzXn03PMHIffv/dpol277apITZm7tWL3ZV1WUXE9yl48YRjx+zLWday/PFH1PEcOeyzhTN4efRR15tOnDiBKVOmmOLbzZs3u46zdsWZFipZsqRHTleCK030zTdAjx52KnTevLbFPzOYIiL+QsGLu3HHENM/bCpHTAM1bmzTQmFhNk0UiWmgxYsXm4Bl0aJFrrRQ8uTJzUwhrrJwxhDvi7ijMTNrwGfMsPefftouCrLUSkTEnyh4ud+XsSyevXDBDj+kSpXsdmaukjBgadnSpoai2bJliwlYuNJy/Phx1/Fy5cqZgKVVq1amP4uIu2zcaH8UObuT8TOnSLzxhtJEIuKfFLwkBIcZcmsz00I7dnCPsq1tofTpgX37bIooGtauONNCv7PgIBI73bZv396khR577LH7+26K3MbZs8BTT9n/FygATJsGPP64LpWI+C8FL/Edq8tRukuXAjdv2mOpUwMMOrjPlP+myMCFTeOYDmLAwvQQm8oRdwc1atTIrLKEhYWZ3UMinpIxI/DZZ8D339tt0bEWAkVE/E7APGtGH8zoMUwPLVli/80KR/Zjad48xlBEjopiwa0zLcT+LE4VKlQwAUvLli2RmdukRTyEu/LZ4p9TJIi1LrxpPpGIBAINZowPBkZ8CcuApVChGG9ip9vJkyeboIUdcJ1y587tSgsVLVrUbd84kTtZtMjWiHNeJzeuaaVFRPyBBjN6Chti9OnjunvlyhV89913JmBhMznnqk/KlCnRpEkTs8pSq1YtJFWPdUlEXBTMnx8IDVVBrogEpoBJGyUWpoXYnp9N5JgW4jRnp0qVKpmApUWLFsiUKZNXz1OCC+vGixSxaSHWuLA8i+OtlCYSkUCk4CWOmBaaOHGiWWXZzul1kfLmzetKCxUuXNhT3yeR22Ld+KBBwLvv2obN3bvb47E2u4mIBBQFL3H0888/oxf7pwNIlSoVnn32WbPKUrNmTaWFxCvYIoithLiLyNnLRUQkGCh4iSNubWa322bNmqF58+bIyLV5ES9ZtQpo3Ro4fJjBNDBkCNCli74dIhIcFLzEUerUqU1Rroi300QDBwIffGD/zToXtvsvUULfFxEJHgpeRPxo3mf79rZHIjFlNGwYkDatt89MRCRxKXgR8QPLlwNt27JwHEiTBhg+3AYvIiLBKIm3T0BE7oytg/r1A2rVsoELJ1GsX6/ARUSCm1ZeRHxYu3Z2kCI9/7zdDs2VFxGRYKaVFxEfxoCFo7MmTwZGjVLgIiJCWnkR8SEcPL5tG1CqlL3/1FPA339rPpGISHRaeRHxEadOATVqAFWrAnv3Rh3XYEURkZgUvIj4CKaHnKIHLyIiEpPSRiJedO0ah30CKVIAyZIBU6cCERFAoUL6toiI+M3Ky/nz51G+fHmULl0aJUqUwChWKYoEoP/9D3jySSByZJaRN68CFxGRewlxOPi6z3fcuHEDV65cQZo0aXDp0iUUL14c69evR5YsWeL0/ufOnTNzh86ePYsM0dfhRXzI/PlA587A6dNApkzAzp2aBC0iwe1cPJ6/fW7lJWnSpCZwoYiICBPM+Fh8JZJgV68CPXoAjRvbwKVCBWDzZgUuIiLxEe/gZdWqVWjQoAFy586NkJAQzJs375bHDB8+HAULFkSqVKlQtmxZrF69Ol6f48yZMyhVqhTy5s2LXr16IWvWrPE9TRGfs28fUKWKbTRHb74J8FfjwQe9fWYiIgEevFy8eNEEFkOHDr3t26dPn44ePXqgb9++2Lx5M6pWrYp69erhwIEDrscwoGE6KPbt8OHD5u2ZMmXC77//jv3792PKlCk4xol0d8AUE5eaot9EfM2sWUCZMsCGDUDmzMCCBcCgQbZQV0REErHmhSsvc+fORWOugUeqWLEiQkNDMWLECNexokWLmscMHDgw3p/jpZdeQs2aNdG8efPbvr1///4YMGDALcdV8yK+gDuHuMLCQYpUubLdUZQ/v7fPTETEt3it5uXq1avYuHEjwsLCYhzn/bVr18bpY3CVxbl6wv8zTVW4cOE7Pr5Pnz7mC3XeDh48eJ9fhYh77NljgxVn4PLOO8DPPytwERHxqT4vJ06cMAW2OXLkiHGc949yJG4c/PPPP+jSpYsp0uXtlVdeQcmSJe/4+JQpU5qbiC9ZvBho2RK4cAFgydaECUC9et4+KxGRwOCRJnVMJ0XHICT2sTthPcyWLVs8cVoiieahh2zzObb6Z5ooTx5dfBERnwxeuCuIW51jr7KEh4ffshrjbsOGDTM3rvyIeMOZM7ZnCxUpYncSlShhO+eKiIj7uLXmJUWKFGblZOnSpTGO835lJv89qHv37ti+fbtpaCeS2CZOBAoUAFaujDrG3UUKXERE3C/erwkvXLiAvdGmxnE7M9M8mTNnRv78+dGzZ0+0b98e5cqVQ6VKlTBy5EizTbpbt27uPncRn7FsGQvMAU6zqFbN22cjIhLY4h28bNiwATVq1HDdZ7BCHTt2xLhx49CyZUucPHkSH374IY4cOWL6tyxevBgF+LJUJEANG8Z6LeDll719JiIigc/nZhu5o+Zl9+7d6vMiHsPfmHHjgO+/B6ZNA5L43JANEZHA7vMSMMGLkwYziidx6/NLLwGTJtn73EnUqpWuuYhIYj5/ax+ESBz98QfARs+7d3OAKPDRR0CLFrp8IiKJTcGLyD1wbXLkSOD11zlLy/ZsYbroiSd06UREvCFgsvWsdylWrBjKly/v7VORAMIdRK1bA9wsx8Clfn2APRQVuIiIeI9qXkTuYNMmmxb66y/br4VzRbm5TgW6IiLup5oXkftME3HrM6dBX71qBykyTVSpki6riIgvUM2LSKwW/126AHPm2PsNGwJjxwKZM+syiYj4ioCpeRFxB47GWrcOSJ4c+PJLYN48BS4iIr4mYFZeNJhREsrZ6YiDz7NkAWbOtFuhVfstIuKbVLArQe3UKaBzZ6BxY/t/ERHx/YJdpY0kqE2YACxYYItzz5/39tmIiEhQpY1EEuK114AdO2wfl/TpdQ1FRPyBVl4kqJw4YQOWy5ftffZs+eYboEwZb5+ZiIjElVZeJGisXm275R46BFy/Dgwf7u0zEhGRoF550XgAuZObN4FPPgGqV7eBS+HCNk0kIiL+SbuNJKAdOwa0bw8sXWrv899ccUmXzttnJiIi0Wk8gAiA5cuBtm2Bo0eB1Klty/9OnWw/FxER8V8BkzYSid4lt39/oFYtG7gUKwZs2GD7uChwERHxfyrYlYBy5IhdbVmxwt5/7jlgyBAgTRpvn5mIiLiLghcJGKxradcOCA8H0qYFvv7a3hcRkcCi4EUCZkfR++/bwKVkSWD6dKBIEW+flYiIeELA1Lxoq3RwY7O5qVOB118Hfv1VgYuISCDTVmnxW4sXA3/8AfTu7e0zERGR+6Wt0hLw/vwTePpp++/KlYEnn/T2GYmISGJRzYv4pRIlgJdfBpImBSpW9PbZiIhIYlLwIn5j4UKgbFkgVy57n1ugWesiIiLBRX/6xeddvQr07Ak0aGB7uLAJHSlwEREJTlp5EZ+2fz/QsiWwfr29X7q03RbNdJGIiAQnBS/is+bMsR1yz54FHngAGDcOaNjQ22clIiLeprSR+JwrV4BXXwWaNrWBS6VKwObNClxERMRS8CI+Ze9eu/V56FB7v1cvYOVKoEABb5+ZiIj4ioAJXtRh1/+xpX9oKLBpE5AlC7BoEfDvfwPJk3v7zERExJeow6543eXLwBtvAN98Y+8/8YRt9Z83r7fPTEREfLHDbsCsvIj/GjDABi4hIUDfvsCKFQpcRETkzhS8iNf16WOLcpcsAT7+GEimPXAiInIXCl4k0V26BHz9NeBw2PsZMwJr1gC1a+ubISIi96bXuJKorl+3NS3c+sxmc5xPREwZiYiIxIVWXiRRMSXUvj2QMydQtKguvoiIxJ+CF/G4CxeAv/6Kut+jB7BtG1Cjhi6+iIjEn4IX8ag//wTKlwfq1wfOn49KEWXOrAsvIiIJo+BFPILFuKNGARUqADt3AhcvAv/7ny62iIjcPwUv4nbnzgFt2gAvvABERNhVly1bgOLFdbFFROT+KXgRt+IuorJlgWnTbHHuZ58B330HZM2qCy0iIu6hrdLitjTRiBG2zf/Vq0D+/DaAYfM5ERERdwqYlRcNZvSeM2eAFi2A7t1t4NKwoV2BUeAiIiKeoMGMcl/WrwdatgT277fTn5kmev11NZ0TERHPDWZU2kgSbO1aoHp14No1oGBBYPp0uy1aRETEkxS8SIJxG3TFikCOHMDo0UCmTLqYIiLieQpeJF42bQIeewxImdLuJlq8GEiXTmkiERFJPAFTsCuex91EXGl5552oY+nTK3AREZHEpeBF4ixvXjsV+tgx4MYNXTgREfEOpY3knt1ynUXfDRoA//d/dvWF84lERES8QSsvcls3bwL/+hdQqBBw4EDU8ccfV+AiIiLepeBFbsG0UN26QN++wPHjwOTJukgiIuI7lDaSGFassEMVjx4FUqdm52KgUyddJBER8R1aeRGDBbgDBgC1atnApVgx2z23c2eliURExLdo5UVw5AjQrh2wfLm9GAxYhgwB0qbVxREREd+j4CXILV1qA5fwcBussJdL+/bePisREZE7U9ooSLFfy3vvAXXq2MClRAlgwwYFLiIi4vu08hKELl0C6tUDVq2y9194AfjyS1ugKyIi4uu08hKE0qSxU6A5k2jqVOCbbxS4iIiI//DZ4OXSpUsoUKAA3nrrLW+fSkC4dg04ezbqPrdAb94MtGrlzbMSEREJoODlk08+QUX2oZf7dvAgUL060Lat7ZxLLM5l91wRERF/45PBy549e7Bz507Ur1/f26cSEE6dAjZuBFav5rX19tmIiIgkcvCyatUqNGjQALlz50ZISAjmzZt3y2OGDx+OggULIlWqVChbtixW81kzHpgqGjhwYHxPTaJxOKL+XaqUrW1hmqhwYV0mEREJsuDl4sWLKFWqFIYOHXrbt0+fPh09evRA3759sXnzZlStWhX16tXDgWjT/RjQFC9e/Jbb4cOHMX/+fDz66KPmJgmzfz9QtSrw229Rx5o0AR56SFdURET8X4jDEf01ejzfOSQEc+fORePGjV3HWKcSGhqKEex2Fqlo0aLmMXFZTenTpw8mTZqEpEmT4sKFC7h27RrefPNNfPDBB7d9/JUrV8zN6dy5c8iXLx/Onj2LDBkyINjMmQM895wtzi1TxqaLQkK8fVYiIiJ3x+fvjBkzxun52601L1evXsXGjRsRFhYW4zjvr127Nk4fgwHOwYMH8ffff2PQoEHo2rXrHQMX5+P5xTpvDFyCEeO3V18Fmja1gcvjjwNz5ypwERGRwOPW4OXEiRO4ceMGcuTIEeM47x/ltD8P4EoNozTnjYFPsNm7F6hcGXBm8nr1sg3oChTw9pmJiIj4SYddppOiY2Yq9rG46NSp0z0fkzJlSnMLVtOnA127AufPA1myABMmANqkJSIigcytKy9Zs2Y1tSqxV1nCw8NvWY1xt2HDhqFYsWIoX748gsHly0C3brbJHAOXJ54AtmxR4CIiIoHPrcFLihQpzE6ipRxVHA3vV2Zew4O6d++O7du3Y/369Qh0u3bZmha29eeC1rvvAitWAHnzevvMREREfDBtxB1Ae1lkEWn//v3YsmULMmfOjPz586Nnz55o3749ypUrh0qVKmHkyJFmm3Q3LhPIfZsyxQ5SvHgRyJYNmDSJBdG6sCIiEjziHbxs2LABNWrUcN1nsEIdO3bEuHHj0LJlS5w8eRIffvghjhw5Yvq3LF682MwpEvf0cGHgwnb/kycDuXPrqoqISHC5rz4vvoQ1L7xxt9Pu3bsDqs8L5xEliUzw3bgBTJwItG8PJE3q7TMTERFJ/D4vARO8JOSL93X8zowfb2tbli8HUqf29hmJiIgEWJM6cS82m+vdG/j1VxvAiIiIiIf6vIh7ZMpk61oYvLz2mq6qiIhIQAUv0Wte/DlNNHq0DVqaN7fHnnrK3kRERMRSzYuPOHcOePFFYNo0IH16YNs2IEjHNImISBA6F4+al4BZefFnmzcDLVrYGUXcQfTee0CePN4+KxEREd+k4MXLaaIRI4A33uBEbrvSwpUXDzcjFhER8WsKXry4k+j554FZs+z9Bg2AsWPtcEUREREJgq3S/jSYkeOXypSxgUuyZMDgwcD8+QpcRERE4kIFu4mcJvrvf4G33wauXQMefBCYPh2oUCExz0JERMT3qGDXB506BTz3nF1hoSZNgG+/tduiRUREJAjTRr5u5UobuKRIAQwZAsyercBFREQkIVSwm0i40vLRR0C9ekDZson1WUVERAKPVl485MQJoGNH4MiRqGPs36LARURE5P4EzMqLr40H6NAB+P57G8QsWuTtsxEREQkc2m3kITt3Am3b2qLcUqU89VlERESCb7eR0kZuEh5uu+M6FSkCbNigwEVERMTdAiZt5E0//wy0aQMcOwbkygVUq2aPh4R4+8xEREQCj1Ze7gPLaz78EHjqKVuYW7gwkDWr+745IiIiciutvCTQ0aO2pmX5cnu/c2fbvyVt2oR+RBEREYkLBS8J8NNPQLt2Nk3EYIWTodu3T8hHEhERkfhS2igerl8H3n8fCAuzgUuJErYoV4GLiIhI4gmY4CUxpkr36AF8/LEdsPjCC8Bvv9ldRSIiIpJ41OclHvbuBapXBwYNAlq18tw3RUREJNici0efF9W8xEOhQsBffwEpU97vt0hEREQQ7GmjxKLARURExLsUvIiIiIhfUfAiIiIifkXBi4iIiPgVBS8iIiLiVxS8iIiIiF8JmOAlMZrUiYiIiPepSZ2IiIj4VZO6gFl5ERERkeCg4EVERET8ioIXERER8SsKXkRERMSvKHgRERERv6LgRURERPxKMgQYh8Ph2nIlIiIi/sH5vO18Hg+q4OX8+fPm//ny5fP2qYiIiEgCnsfZ7yWomtTdvHkThw8fRvr06RESEuKxz8NOvuvXr4evSIzzcefnuN+PldD3j+/7xfXx93ocX1EwoD548OA9my8FEv2eePf66ffEP/jS70n5RDqX230ehiMMXHLnzo0kSZIE18oLv+C8efN6/PMkTZrUp56EEuN83Pk57vdjJfT94/t+cX18XB/Hx/jSz42n6ffEu9dPvyf+wZd+T5Im0rnc6fPca8XFSQW7CdS9e3cE2/m483Pc78dK6PvH9/3i+nhf+3nwFb52XfR74pnrpN+TwPk96Z5I53K/nyfg0kYi/j6zQyRY6fdE4korLyKJIGXKlOjXr5/5v4jo90Tuj1ZeRERExK9o5UVERET8ioIXERER8SsKXkRERMSvKHgRERERv6LgRURERPyKghcRH8MRAtWrV0exYsVQsmRJzJw509unJOKTmjRpggceeADNmjXz9qlIItNWaREfc+TIERw7dgylS5dGeHg4QkNDsWvXLqRNm9bbpybiU1asWIELFy5g/PjxmDVrlrdPRxKRVl5EfEyuXLlM4ELZs2dH5syZcerUKW+flojPqVGjhhnCK8FHwYtIPK1atQoNGjQwk085uXzevHm3PGb48OEoWLAgUqVKhbJly2L16tUJus4bNmwwk9I5kVrEnyTm74kEHwUvIvF08eJFlCpVCkOHDr3t26dPn44ePXqgb9++2Lx5M6pWrYp69erhwIEDrsfwD3Xx4sVvuR0+fNj1mJMnT6JDhw4YOXKkvkfidxLr90SCk2peRO7nFygkBHPnzkXjxo1dxypWrGjqVEaMGOE6VrRoUfOYgQMHxunjXrlyBbVr10bXrl3Rvn17fY/Er3nq94R+/vlnEyCp5iW4aOVFxI2uXr2KjRs3IiwsLMZx3l+7dm2cPgYHvXfq1Ak1a9ZU4CIByR2/JxLcFLyIuNGJEydw48YN5MiRI8Zx3j969GicPsaaNWvMkjprBFi4y9uff/6p75MEDHf8nlCdOnXQvHlzLF68GHnz5sX69es9cLbii5J5+wREAnWZPPZqSuxjd/LEE0+YIl2RQHc/vye0ZMkSD5yV+AOtvIi4UdasWZE0adJbXj2yX0vsV5kiwUq/J3K/FLyIuFGKFCnMDomlS5fGOM77lStX1rUW0e+JuIHSRiLxxI6ee/fudd3fv38/tmzZYprJ5c+fHz179jSFtuXKlUOlSpXMVmdu/+zWrZuutQQN/Z6IRzlEJF5WrFjh4K9O7FvHjh1djxk2bJijQIECjhQpUjhCQ0MdK1eu1FWWoKLfE/Ek9XkRERERv6KaFxEREfErCl5ERETEryh4EREREb+i4EVERET8ioIXERER8SsKXkRERMSvKHgRERERv6LgRURERPyKghcRERHxKwpeRMRvrF27FiEhIahbt663T0VEvEjjAUTEbzz//PO4dOkSZs+ejT179phBmCISfLTyIiJ+4eLFi5g+fTp69OiBmjVrYty4cd4+JRHxEgUvIuIXGLjkzJkTFSpUQNu2bTF27Fg4HBzoLSLBRsGLiPiFMWPGmKCFGjdujPDwcCxbtszbpyUiXqCaFxHxebt27UKRIkXM/x999FFzrE2bNub/U6ZM8fLZiUhi08qLiPjFqkv58uVdgQtxFWbOnDk4ffq0V89NRBKfghcR8WnXr1/HhAkTXCstTnXq1EH69OkxefJkr52biHhHMi99XhGROFm4cCGOHTuG4sWLY+vWrTHeVrVqVbMq88orr+hqigQR1byIiE9r0KCBCWDuZuPGjQgNDU20cxIR71LwIiIiIn5FNS8iIiLiVxS8iIiIiF9R8CIiIiJ+RcGLiIiI+BUFLyIiIuJXFLyIiIiIX1HwIiIiIn5FwYuIiIj4FQUvIiIi4lcUvIiIiIhfUfAiIiIifkXBi4iIiMCf/D9MN9MLFNo/+QAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def right_sided_sum(f, a, b, n):\n", " \"\"\"Approximate the integral of f over [a, b] using \n", " the right-sided sum and n evenly spaced points.\"\"\"\n", " x = np.linspace(a, b, n+1)[1:]\n", " return np.sum(f(x))*(x[1]-x[0])\n", "\n", "# reference integral value\n", "I_ref = np.exp(1) - 1.0 \n", "# number of points to use in each sum\n", "n_vals = 2**np.arange(1, 10)\n", "# initialize the approximate integrals\n", "I_appx = np.zeros(len(n_vals))\n", "for i in range(len(n_vals)):\n", " I_appx[i] = right_sided_sum(np.exp, 0, 1, n_vals[i])\n", "# absolute error\n", "I_err = abs(I_appx-I_ref)\n", "# Delta = (b-a)/(n+1)\n", "Delta = 1/(n_vals+1)\n", "# Plot the error along with linear and quadratic forms\n", "plt.loglog(Delta, I_err, 'k', label='error')\n", "plt.loglog(Delta, Delta, 'r--', label='$\\mathcal{O}(\\Delta)$')\n", "plt.loglog(Delta, 10*Delta**2, 'b-.', label='$\\mathcal{O}(\\Delta^2)$')\n", "plt.xlabel('$\\Delta$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Surely, the error is much closer to the $\\mathcal{O}(\\Delta)$ term and, hence, the right-sided sum is *first order*.\n", "\n", "
\n", "\n", "**Exercise**: Determine the order of the left-sided sum, midpoint rule, and trapezoid rule. You may use the definite integral $\\int^1_0 e^x dx = e - 1$ as a test case. \n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using `quad` for Callable Functions\n", "\n", "The previous section started with samples. A different problem is: \"I have a function I can evaluate anywhere, and I want the integral.\" That is when `quad` is usually the first SciPy tool to try.\n", "\n", "SciPy is a large collection of open-source, Python tools for numerical and other scientific computing. For example, NumPy is a major component of SciPy but is often used as a standalone package (as are Matplotlib and SymPy). The `scipy.integrate` module contains several integration routines, but `quad` is the workhorse for one-dimensional integrals of callable functions." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2026-07-02T15:33:04.153787Z", "iopub.status.busy": "2026-07-02T15:33:04.153682Z", "iopub.status.idle": "2026-07-02T15:33:04.155968Z", "shell.execute_reply": "2026-07-02T15:33:04.155591Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "quad result = 1.7182818284590453\n", "estimated error = 1.91e-14\n", "true error = 2.22e-16\n" ] } ], "source": [ "from scipy.integrate import quad\n", "\n", "I_quad, err_quad = quad(np.exp, 0.0, 1.0)\n", "I_ref = np.exp(1.0) - 1.0\n", "\n", "print(f\"quad result = {I_quad:.16f}\")\n", "print(f\"estimated error = {err_quad:.2e}\")\n", "print(f\"true error = {abs(I_quad - I_ref):.2e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, `quad` returns both the approximate integral and an estimated absolute error. Here, that error is ridiculously small (close to machine precision). For many smooth one-dimensional problems, that's all one needs to do in order to compute an integral.\n", "\n", "That does not mean `quad` is magic. Under the hood, `quad` calls routines from QUADPACK. For ordinary finite intervals, the method adaptively subdivides the interval and applies Gauss-Kronrod-style quadrature rules on the subintervals. It keeps subdividing until the estimated error satisfies tolerances such as `epsabs` and `epsrel`, or until it runs into a subdivision limit.\n", "\n", "A few practical rules help:\n", "\n", "- Use `quad` for callable functions, not for data that are already sampled.\n", "- If the integrand has a jump, cusp, singularity, or known trouble spot, split the integral or use the `points` argument.\n", "- Infinite intervals are allowed; internally, they are transformed to finite intervals.\n", "- The returned error is an estimate, not a proof. A physical or graphical sanity check is still your job." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2026-07-02T15:33:04.156912Z", "iopub.status.busy": "2026-07-02T15:33:04.156820Z", "iopub.status.idle": "2026-07-02T15:33:04.159180Z", "shell.execute_reply": "2026-07-02T15:33:04.158651Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.09998291963330941 1.1161971513740099e-08\n", "0.09998291965798728 1.15733619018813e-13\n" ] } ], "source": [ "# A known internal point can help quad spend effort in the right place.\n", "# This toy function is continuous, but it has a sharp change in slope at x=0.4.\n", "def peaked(x):\n", " return np.exp(-20*abs(x - 0.4))\n", "\n", "I_plain, err_plain = quad(peaked, 0.0, 1.0)\n", "I_hint, err_hint = quad(peaked, 0.0, 1.0, points=[0.4])\n", "\n", "print(I_plain, err_plain)\n", "print(I_hint, err_hint)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiple Integrals and Vector-Valued Integrands\n", "\n", "Other SciPy integration functions are useful when the geometry gets more interesting:\n", "\n", "- `dblquad` computes a double integral with inner limits that may depend on the outer variable.\n", "- `nquad` generalizes the same idea to more variables.\n", "- `quad_vec` is useful when the integrand returns a vector rather than a scalar.\n", "\n", "The argument order for `dblquad` is easy to trip over: the integrand is written as `f(y, x)`, with the inner variable first." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2026-07-02T15:33:04.160249Z", "iopub.status.busy": "2026-07-02T15:33:04.160155Z", "iopub.status.idle": "2026-07-02T15:33:04.163509Z", "shell.execute_reply": "2026-07-02T15:33:04.163103Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "J_in = -9.42477796\n", "-3*pi = -9.42477796\n", "estimated error = 6.00e-09\n" ] } ], "source": [ "from scipy.integrate import dblquad\n", "\n", "# Current balance example: if J(x,y) = (x, 2y), then div J = 3.\n", "# The net current into the unit disk is -int_D div(J) dA = -3*pi.\n", "def minus_divergence(y, x):\n", " return -3.0\n", "\n", "J_in, err = dblquad(\n", " minus_divergence,\n", " -1.0, 1.0,\n", " lambda x: -np.sqrt(1.0 - x**2),\n", " lambda x: np.sqrt(1.0 - x**2),\n", ")\n", "\n", "print(f\"J_in = {J_in:.8f}\")\n", "print(f\"-3*pi = {-3*np.pi:.8f}\")\n", "print(f\"estimated error = {err:.2e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the same accounting idea that shows up in transport: local rates are useful, but many engineering answers are totals over a region, a surface, an energy interval, or a time interval." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical Probability Accounting\n", "\n", "Numerical integration is often the step between a local model and a quantity we can check. The example below is deliberately self-contained: a particle starts uniformly in a cell of width $h$, and an event occurs along distance with rate $\\beta$. If the particle starts at position $x$, the probability of reaching the left face before the event is $e^{-\\beta x}$, and the probability of reaching the right face before the event is $e^{-\\beta(h-x)}$.\n", "\n", "For an isotropic initial direction, the two face probabilities are\n", "\n", "$$\n", "P_L = \\frac{1}{h}\\int_0^h \\frac{1}{2}e^{-\\beta x}\\,dx,\n", "\\qquad\n", "P_R = \\frac{1}{h}\\int_0^h \\frac{1}{2}e^{-\\beta(h-x)}\\,dx.\n", "$$\n", "\n", "The complementary event probability is\n", "\n", "$$\n", "P_C = 1-P_L-P_R.\n", "$$\n", "\n", "A good numerical habit is to carry the denominator in words even when the result is dimensionless:\n", "\n", "```text\n", "P_L [left-face crossings / particle born in the cell]\n", "P_R [right-face crossings / particle born in the cell]\n", "P_C [events before face crossing / particle born in the cell]\n", "```\n", "\n", "The accounting check is $P_L+P_R+P_C=1$.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import quad\n", "\n", "\n", "def face_probabilities_quad(beta, h):\n", " \"\"\"Return left, right, and collision/event probabilities by quadrature.\"\"\"\n", " left_integrand = lambda x: 0.5*np.exp(-beta*x)\n", " right_integrand = lambda x: 0.5*np.exp(-beta*(h - x))\n", "\n", " P_left, err_left = quad(left_integrand, 0.0, h)\n", " P_right, err_right = quad(right_integrand, 0.0, h)\n", "\n", " P_left /= h\n", " P_right /= h\n", " P_event = 1.0 - P_left - P_right\n", " err = (err_left + err_right)/h\n", " return P_left, P_right, P_event, err\n", "\n", "\n", "def face_probability_closed(beta, h):\n", " \"\"\"Closed-form one-face probability for beta*h not too small.\"\"\"\n", " tau = beta*h\n", " if tau == 0.0:\n", " return 0.5\n", " return 0.5*(1.0 - np.exp(-tau))/tau\n", "\n", "\n", "for tau in [0.01, 0.1, 1.0, 10.0]:\n", " beta = tau\n", " h = 1.0\n", " P_left, P_right, P_event, err = face_probabilities_quad(beta, h)\n", " P_face = face_probability_closed(beta, h)\n", " print(\n", " f\"tau={tau:5.2f} \"\n", " f\"P_left={P_left:.6f} P_right={P_right:.6f} \"\n", " f\"P_event={P_event:.6f} sum={P_left + P_right + P_event:.12f} \"\n", " f\"closed_face={P_face:.6f}\"\n", " )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same integral can be evaluated with a fixed quadrature rule. This is useful when you want full control over the points and weights or when the same rule will be reused many times.\n", "\n", "For Gauss-Legendre quadrature, NumPy gives points and weights on $[-1,1]$. To integrate over $[0,1]$, map\n", "\n", "$$\n", "\\eta_q = \\frac{z_q+1}{2},\n", "\\qquad\n", "w_q^{[0,1]} = \\frac{w_q^{[-1,1]}}{2}.\n", "$$\n", "\n", "Then\n", "\n", "$$\n", "\\int_0^1 g(\\eta)\\,d\\eta \\approx \\sum_q w_q g(\\eta_q).\n", "$$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from numpy.polynomial.legendre import leggauss\n", "\n", "\n", "def unit_gauss_legendre(n):\n", " \"\"\"Gauss-Legendre points and weights mapped to [0, 1].\"\"\"\n", " z, w = leggauss(n)\n", " eta = 0.5*(z + 1.0)\n", " weights = 0.5*w\n", " return eta, weights\n", "\n", "\n", "def face_probability_fixed_quad(beta, h, n):\n", " \"\"\"One-face probability using an n-point rule on eta in [0, 1].\"\"\"\n", " eta, weights = unit_gauss_legendre(n)\n", " values = 0.5*np.exp(-beta*h*eta)\n", " return np.sum(weights*values)\n", "\n", "\n", "beta = 3.0\n", "h = 0.7\n", "exact = face_probability_closed(beta, h)\n", "\n", "print(f\"closed form = {exact:.14f}\")\n", "for n in [1, 2, 3, 4, 8, 16]:\n", " approx = face_probability_fixed_quad(beta, h, n)\n", " print(f\"n={n:2d} approx={approx:.14f} abs_error={abs(approx - exact):.2e}\")\n", "\n", "eta, weights = unit_gauss_legendre(4)\n", "print(f\"weight sum check = {weights.sum():.16f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further Reading\n", "\n", "- The SciPy documentation for `scipy.integrate` is the best practical reference for the current API.\n", "- The `quad` documentation is worth reading once, especially the notes about QUADPACK routines, tolerances, `points`, and infinite intervals.\n", "- The sampled-data tools to remember are `trapezoid`, `simpson`, and `cumulative_trapezoid`. Older names such as `trapz`, `simps`, and `cumtrapz` should not be used in new SciPy code." ] } ], "metadata": { "celltoolbar": "Edit Metadata", "kernelspec": { "display_name": "Python 3", "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.13" } }, "nbformat": 4, "nbformat_minor": 2 }