# 积分 ## 符号积分 积分与求导的关系: $$\frac{d}{dx} F(x) = f(x) \Rightarrow F(x) = \int f(x) dx$$ 符号运算可以用 `sympy` 模块完成。 先导入 `init_printing` 模块方便其显示: In [1]: ```py from sympy import init_printing init_printing() ``` In [2]: ```py from sympy import symbols, integrate import sympy ``` 产生 x 和 y 两个符号变量,并进行运算: In [3]: ```py x, y = symbols('x y') sympy.sqrt(x ** 2 + y ** 2) ``` Out[3]:$$\sqrt{x^{2} + y^{2}}$$ 对于生成的符号变量 `z`,我们将其中的 `x` 利用 `subs` 方法替换为 `3`: In [4]: ```py z = sympy.sqrt(x ** 2 + y ** 2) z.subs(x, 3) ``` Out[4]:$$\sqrt{y^{2} + 9}$$ 再替换 `y`: In [5]: ```py z.subs(x, 3).subs(y, 4) ``` Out[5]:$$5$$ 还可以从 `sympy.abc` 中导入现成的符号变量: In [6]: ```py from sympy.abc import theta y = sympy.sin(theta) ** 2 y ``` Out[6]:$$\sin^{2}{\left (\theta \right )}$$ 对 y 进行积分: In [7]: ```py Y = integrate(y) Y ``` Out[7]:$$\frac{\theta}{2} - \frac{1}{2} \sin{\left (\theta \right )} \cos{\left (\theta \right )}$$ 计算 $Y(\pi) - Y(0)$: In [8]: ```py import numpy as np np.set_printoptions(precision=3) Y.subs(theta, np.pi) - Y.subs(theta, 0) ``` Out[8]:$$1.5707963267949$$ 计算 $\int_0^\pi y d\theta$ : In [9]: ```py integrate(y, (theta, 0, sympy.pi)) ``` Out[9]:$$\frac{\pi}{2}$$ 显示的是字符表达式,查看具体数值可以使用 `evalf()` 方法,或者传入 `numpy.pi`,而不是 `sympy.pi` : In [10]: ```py integrate(y, (theta, 0, sympy.pi)).evalf() ``` Out[10]:$$1.5707963267949$$In [11]: ```py integrate(y, (theta, 0, np.pi)) ``` Out[11]:$$1.5707963267949$$ 根据牛顿莱布尼兹公式,这两个数值应该相等。 产生不定积分对象: In [12]: ```py Y_indef = sympy.Integral(y) Y_indef ``` Out[12]:$$\int \sin^{2}{\left (\theta \right )}\, d\theta$$In [13]: ```py print type(Y_indef) ``` ```py ``` 定积分: In [14]: ```py Y_def = sympy.Integral(y, (theta, 0, sympy.pi)) Y_def ``` Out[14]:$$\int_{0}^{\pi} \sin^{2}{\left (\theta \right )}\, d\theta$$ 产生函数 $Y(x) = \int_0^x sin^2(\theta) d\theta$,并将其向量化: In [15]: ```py Y_raw = lambda x: integrate(y, (theta, 0, x)) Y = np.vectorize(Y_raw) ``` In [16]: ```py %matplotlib inline import matplotlib.pyplot as plt x = np.linspace(0, 2 * np.pi) p = plt.plot(x, Y(x)) t = plt.title(r'$Y(x) = \int_0^x sin^2(\theta) d\theta$') ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXIAAAEVCAYAAAD91W7rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4U+W1x/HvYlBUFFQUFVDqgMOtKFARxdagSAEtVkUFZ7FXbatWqx2w3Epbr1VbxaF1QrAgjq0IDogWL1FBZZBBEFGEqjgAUsSCiAJn3T/egDEcTpJzkrOzk9/nefKcDG92Fvs5rLxn7Xcwd0dEROKrQdQBiIhI3SiRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEXKOoAxApBjNrBiSBVcBxwM7AgcCxwFJgnrv/M7IARQpIPXIpV78CngeaA/sCuPvzhIR+F/B/0YUmUlhK5FJ2zKwBcB7wN+B77v4GsNrMWgLLgK2B7SILUKTAlMilHB0ObHD3ue6+MvXcICABfA4c4+7/iSo4kUJTjVzKUTdgcvoT7v7LiGIRKTr1yKUcHQ1MizoIkfqiRC5lxcwaAkcAM6OORaS+KJFLQZnZt3Jos7uZbVukEA4FmgKzi3T8GpnZfmZ2kpldY2Yds7Td4rkq8jmSMqNELgVjZnsDXXJo+glQrJr1EcBSd/93kY6fzQnAh8DNwFVbapTDuSrmOZIyo0QuhXSxuz+UrZG7rweeNrNzihDD4cDcIhw3J+4+xN2nAm2Af9XQ9KKazlWRz5GUGSVyyZmZfcfMnjCzj82sZ9rzB5vZp8AFuZYD3H0a0L0IYXahCInczK43sx55vOUk4H+3cKxDgA8ynutqZgPM7Coz2xFqPkdm1s3MPjKzNjUdQyqDErnkzN2nA5cCOwJT0176L+DPhMk3a/I45Cdmtm+h4jOznYB9gDmFOuZG7v5rd38uxzj6ALcBrbbQ5ATSZpamzsF57j4ceA/om9Z2S+foJUIJaXEOx5Ayp3Hkkhd3f8/MXgHOBG43s17A20A/4I95Hm420Al4Z+MTqdrxf9fwnlfdfewWXuuU+lnwRJ4rMzsJuJrwhZek+l75YcB1aY9vSHu8P7Ah7bXNzlFKJ745xLKmY0iZUyKX2hgBXGJm04Gv3H2GmW3r7lXpjVI90w3AdwnJtSfwv+4+P9XkU6Bd+nvcfREwsJZxdQKqqENpxcx2AM4BFgEHEKb59wBOdvfTzKwToQ6/BzAdaAgc7+4DUvE/Djye5WO2dXdPfd4ehMTeycy+Q/iC/ENa203nyMzaAWcTJjtdDDyU4zGkzKm0IrXxGHAQ8K3UQlQQEtomZrYnYYXBpwmrDz4NPAK8n9bsC2CrAsbVAXjb3b+owzF+SLhQ+TLh39gBGA9sHCq4CzAf+C93HwOMJkxAykf6uToGeMrd7wFGEcox6SWcL4CtzGw74FHgJncfD+zA1+WZbMeQMqceudRGY2CNuz+Y9tz69Abu/j5AaqGqVak1T57KOE4zYEX6E3UsrRxC3Wd0PkNIzHMIfz08b2aXE3rmuPt4M/sjcH+q/RG1+Mz0c9UaeDN1vw/wtLsvT3t94zk6GZjj7ivNrAnQ1N0/yfEYUuaUyKU2jgYmZTy3xMyauvtqADM7gLDKYEfgxdRzJ7h7ejLfna8TEFD70oqZbUNYrvaufN+bdozDgQvc/QIz25VwwfIuoD/Qw8yOT/2F0Q24PvW2c4ChZtYz1VPORfq5+iR8tBmhbHJxRtvdCX8BHMzXk5yOA15N+8xsx5Ayp9KK5MXMTgcuAxqbWee0l14A0h/3IIzOMKBJ6iLgsozDHUrG4lZ1cCDh97kuPfJlwGup2v4ZwJWp5xcR/i1TUsMrV7r7Z6nXPgd2JeMviyzSz9XfgfbAj4Dfb/xLJs2hhC/Nh4DWqYvLLYB1wPY5HkPKnKWuuYjUiZk1B65y90E5tm8CXOfuPy/Q558NDAd2qGONvOhyPVeFPkdSvmrskZtZEzObYmazzGxeqjaY2SZhZp+Z2czULaf/yFJeUjXw5WbWIse39APuLmAIBwKzSz2JQ17nqtDnSMpUjYnc3dcC3dz9UMKfbt3M7Khqmr7g7h1St2uLEajEwq2EGY01Ss1G/NTd3yrgZx9MmCQTFzWeqyKdIylTWS92ps3U24owbKq6WqAVMiiJp9TY6KE5tFsMLC7wxx9MGN8eC9nOVZHOkZSprBc7zayBmc0i7Dw+0d3nZTRx4Egzm21m48zsoGIEKrIlqan5e1K4C6cisZI1kbt7Vaq00hr4npklMprMANq4+yHA7cCYgkcpUrPvAG+5+8dRByIShbxGrZjZ/wBfuPufa2jzL6CTu2dO9NDwGBGRWnD3GsvX2UattEgNldo44eI4MrbQMrOWqYkIpMYVW2YSTwsmtrdrrrkm8hgqMfZc4j/ggAOYOHFi5HFW6vkv9Vvc489FtouduwMjzKxBKunf72HK8kWpxHw3YbnMH5vZemANYciUSFFNnDiRe+65h+7du9O0aVMSicSm11asWMHQoUPZddddad++PZ06ddrygUTKQI2J3N3nEKZYZz5/d9r9vwJ/LXxoIlvWrl07pk+fzpw5cxg9evQ3XhsxYgTdunWjY8eOnHvuuTzwwAMRRSlSP7TWSo7Se3xxE+fYofr4W7VqxYIFC6ptv2jRIvr27UujRo1YsSKfmfPFUY7nP07iHn8utNZKjuL8yxDn2CH/+KuqqmjYMKwUm7p8E6lKO/+lJu7x50KJXMrO/vvvz9KlS1m7di077LBD1OGIFF29LZplZl5fnyWV7d///jfDhw+nWbNmHHzwwRxxxBFRhyRSa2aGZxl+qEQuIlLCcknkKq2IiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjMaWMJEalY7lBV9fVt42OAbbeNNrZ8KJGLSNlwh+XLYeFCWLQo3BYuhCVLYNWqzW9ffgkNGoSb2df3DzgAZsyI+l+TOy1jKyKx9fnnMHkyTJwIySTMnQtbbw177/3N2x57wPbbb35r0iQk8FKm9chFpKy4h8Q9fnxI3rNnQ8eOkEhAt27QoQM0bx51lIWlRC4iZWH5chg5Eu65J5Q+Tj45JO8jj4xXLbs2cknkqpGLSElyD+WSoUNh3Dg48US4917o2rX0yyH1rcYeuZk1AV4Atga2Asa6+8Bq2t0G9ALWAOe5+8xq2qhHLiJZucOTT8Kvfw0NG8KFF8JZZ8GOO0YdWTTq3CN397Vm1s3d15hZI2CSmR3l7pPSPqQ3sK+772dmhwN3Al0K8Q8QkcoycyZceSUsXQo33QQ9e6r3nYusE4LcfU3q7lZAQ2BFRpM+wIhU2ylAczNrWcggRaS8ffghnH8+9OoFp50WLmL26qUknqusidzMGpjZLGApMNHd52U0aQUsTnv8AdC6cCGKSLn66iv43e+gfXvYbTd4+224+GJopKt3ecl6uty9CjjUzJoBz5pZwt2TGc0yvzerLYYPHjx40/1EIkEikcgnVhEpIwsXQv/+sMsu8Npr0LZt1BGVhmQySTKZzOs9eQ0/NLP/Ab5w9z+nPXcXkHT3h1OP5wNHu/vSjPfqYqeIAPDgg/Czn8GgQXDZZSqh1KTOFzvNrAWw3t1Xmtk2wHHA7zKaPQFcAjxsZl2AlZlJXEQEYPVquPRSePllePbZMJlH6i5baWV3YISZNSDU0+939+fN7CIAd7/b3ceZWW8zewf4HDi/uCGLSBzNmgX9+kGXLqGU0rRp1BGVD83sFJGie/RR+OlP4ZZb4Mwzo44mXjSzU0Qid9ttcOONMGECHHJI1NGUJyVyESkKdxg4EMaMgUmTNCqlmJTIRaTg1q2DCy6ABQvCaoU77xx1ROVNiVxECmrVKujbN6wL/vzz5b86YSnQnp0iUjDLl4d1wffcE0aPVhKvL0rkIlIQK1dCjx7QvXtYN1zT7OuPhh+KSJ2tXh2SeOfOMGSIZmoWknYIEpGi++ILOP542Gef0BNXEi8sJXIRKaqvvoKTTgr7ZI4cGTaCkMJSIheRolm/Pky537AhzNxs3DjqiMqTZnaKSFFUVcGAAaE2PnasknjUlMhFJG9XXQXvvQfPPBPGi0u0lMhFJC933hkS+Msva5x4qVCNXERy9uyzcN55Ye2UffaJOprKoBq5iBTM3Llw9tnw+ONK4qVGMztFJKslS+CEE8J64l27Rh2NZFIiF5EarVkDJ54I558PZ5wRdTRSHdXIRWSLqqrg9NOhSZMw4UezNuufauQiUifXXAMffxyWo1USL11K5CJSrTFjYMQImD5dY8VLnUorIrKZt96C734XnnoqrGgo0cmltKKLnSLyDatWhYWwrrtOSTwu1CMXkU3c4dRTYaedwpK0Er0698jNrI2ZTTSzN8xsrpldVk2bhJl9ZmYzU7dBdQ1cRKJx442weDHcfnvUkUg+sl3sXAdc4e6zzKwp8JqZ/dPd38xo94K79ylOiCJSH/75T7j1Vpg6VRc346bGHrm7L3H3Wan7q4E3gT2qaaqBSSIx9u67Yfr9gw9C69ZRRyP5yvlip5m1BToAUzJecuBIM5ttZuPM7KDChScixfbll9C3L/zyl5BIRB2N1EZO48hTZZV/AD9L9czTzQDauPsaM+sFjAHaVXecwYMHb7qfSCRI6LdGJHJXXgl77QVXXBF1JAKQTCZJJpN5vSfrqBUzaww8BTzj7rdkPaDZv4BO7r4i43mNWhEpMY88Ar/5Dbz2GjRrFnU0Up06T9E3MwOGAfO2lMTNrCWwzN3dzDoTvhxWVNdWRErH22/DJZeENcaVxOMtW2mlK3AW8LqZzUw9dzWwJ4C73w30BX5sZuuBNUC/IsUqIgXyxRehLn7ttdCxY9TRSF1pQpBIBfrRj0IyHzVKi2GVOq1+KCKbGTECJk+GadOUxMuFeuQiFWTuXOjWDSZOhG9/O+poJBdaNEtENlm9Oqyj8qc/KYmXG/XIRSqAO5xzDjRuDMOHRx2N5EM1chEBQvKeOTOsoyLlRz1ykTL3+utw7LHw4otw4IFRRyP5Uo1cpMKtWhXq4jfdpCReztQjFylT7nDWWdCkCQwbFnU0UluqkYtUsHvvDWWVKZnrlUrZUY9cpAzNng3du8NLL8EBB0QdjdSFauQiFeg//wl18SFDlMQrhXrkImVk4+bJLVrAXXdFHY0UgmrkIhVmyJCwbduoUVFHIvVJPXKRMjFpEpxySri42bZt1NFIoahGLlIhli6Ffv3gvvuUxCuRErlIzK1fH5L4+edD795RRyNRUGlFJOYGDoTp02H8eGjYMOpopNB0sVOkzD3xBDzwQNg8WUm8cimRi8TUggVhy7axY2GXXaKORqKkGrlIDH32GfTpA3/4AxxxRNTRSNRUIxeJmQ0b4MQTYc894Y47oo5Gik3DD0XK0KBB8PnncOutUUcipUI1cpEYeegheOSRsNNP48ZRRyOlosYeuZm1MbOJZvaGmc01s8u20O42M1tgZrPNrENxQhWpbNOnw2WXwZgxYS0VkY2y9cjXAVe4+ywzawq8Zmb/dPc3NzYws97Avu6+n5kdDtwJdCleyCKVZ8kSOPlkuPtuaN8+6mik1NTYI3f3Je4+K3V/NfAmsEdGsz7AiFSbKUBzM2tZhFhFKtKXX4YkPmBA+CmSKeeLnWbWFugAZO430gpYnPb4A6B1XQMTEaiqgrPPhtat4be/jToaKVU5XexMlVX+Afws1TPfrEnG42rHGQ4ePHjT/UQiQSKRyClIkUrkDpdfDsuWhen3DTTGrCIkk0mSyWRe78k6jtzMGgNPAc+4+y3VvH4XkHT3h1OP5wNHu/vSjHYaRy6ShxtuCNPvX3wRmjePOhqJSp3HkZuZAcOAedUl8ZQngHNS7bsAKzOTuIjkZ8QIuPNOeOYZJXHJrsYeuZkdBbwIvM7X5ZKrgT0B3P3uVLu/AD2Bz4Hz3X1GNcdSj1wkB+PHw3nnwcSJcOCBUUcjUculR64p+iIlZNq0sKb42LFw5JFRRyOlQFP0RWLkrbfCQljDhimJS36UyEVKwBtvwDHHwPXXh2Qukg+ttSISsdmzoWdP+POf4cwzo45G4kiJXCRCr70Gxx8Pt98Op54adTQSV0rkIhF59dVQRrnnHvjhD6OORuJMiVwkApMmhXVT7rsv9MhF6kIXO0Xq2XPPwUknwahRSuJSGErkIvXEHYYMgXPPhdGjoUePqCOScqHSikg9WLsWLr4YZs2CV16Btm2jjkjKiXrkIkX28ceQSMCaNTB5spK4FJ4SuUgRTZ0Khx0GJ5wQ9trcbruoI5JypNKKSBG4w9ChYcf7oUPhxBOjjkjKmRK5SIG98w5ceCGsWgXJJBx0UNQRSblTaUWkQNavD9Psu3QJwwpfeUVJXOqHeuQiBTB7NlxwATRrBlOmwD77RB2RVBL1yEXqYOlS+MUv4Ljj4Mc/hgkTlMSl/imRi9TC++/DpZeGHXzWrAnjwy+4AKzG5f9FikOJXCQPb70FAwZAhw6w7bYwbx789a+wxx5RRyaVTDVykSxWrYJx48I48EmT4JJLwsiUHXeMOjKRQIlcpBorVsATT4Q1UZJJ6NoVTjkFRo6Epk2jjk7km7T5slS8DRtCD3vmzFDrnjo1bPhw7LEheR9/PDRvHnWUUqly2XxZiVzKWlUVfP45fPJJGGGybFm4LV0KH3wAr78ebrvuGurehx4KHTuGtVE0nV5KgRK5lLSqKli4EBYv/maC3Xh/9Wr46qvNb+vXh/dv/HXa+HPDhs3bbtgQLkruuuvXt5Ytw8/dd4f27eGQQ9TjltJVkERuZsOB44Fl7n5wNa8ngLHAotRTj7n7tdW0UyKvYGvXwty5oXSxsYTx+uvQokVYDTAzye6yC+ywA2y11ea3hg2/HuaX/rNBA9h66y23FYmjQiXy7wKrgZE1JPKfu3ufLMdRIq8w6RcMJ06Evff+unzRoYN6wiK5yCWRZx214u4vmVnbbJ+VR1xSxpYsgTFj4LHHwlT17t3htNPCaA8lbZHiKMTwQweONLPZwIfAVe4+rwDHlRh54w344x/h6aehd++wG86YMbpgKFIfCpHIZwBt3H2NmfUCxgDtqms4ePDgTfcTiQSJRKIAHy9RmjYNrrsurPR3+eVhlmOzZlFHJRJfyWSSZDKZ13tyGrWSKq08WV2NvJq2/wI6ufuKjOdVIy8jyWRI4PPnh0WjLrggjA4RkcIqSI08hw9pSRjR4mbWmfDlsCLb+ySeliwJU9RnzoTf/AbOOiuMDhGR6GRdNMvMHgJeBvY3s8VmNsDMLjKzi1JN+gJzzGwWcAvQr3jhSlTcYfjwMO66XbtQEx8wQElcpBRoQpBktWhR2Lrs009h2LAwfFBE6kcupRUtYytbtGED3HwzdO4M3/9+GE6oJC5SerT6oVTr00/hjDPCpgmvvgr77ht1RCKyJeqRy2bmz4fDD4cDDoDnn1cSFyl1SuTyDU8/Dd/7HgwcCEOGQCP9zSZS8vTfVIAwKuWGG+D222HsWDjiiKgjEpFcKZELa9aECT0LF4ZNFVq1ijoiEcmHSisVbtUq6NkzLAH7wgtK4iJxpERewT77DHr0gIMOgvvvh222iToiEakNJfIKtWJFWGL2sMPgzjtDj1xE4kn/fSvQ8uVhY+FEAm69VTvoiMSdEnmFWbo0JPDeveHGG5XERcqBEnkF+eijkMRPOw2uvVZJXKRcaNGsCvHpp3DUUXDmmXD11VFHIyK5KsjmywUMRok8Il98EUandO4MN90UdTQikg8lcmHDBujbNwwtHDVKo1NE4qZedgiS0uUOP/0prF4NjzyiJC5SrpTIy9i114Yp98mkdvIRKWdK5GXq3nvhb3+DyZNhhx2ijkZEikk18jL05JNha7YXX4T99os6GhGpC13srEBz5sAxx4R1xTt3jjoaEakr7dlZYZYvhxNPDNPulcRFKod65GVi3bowVvzww+H666OORkQKRaWVCnLJJfDuu2F3n4YNo45GRAqlIKUVMxtuZkvNbE4NbW4zswVmNtvMOtQmWKm9oUPDJskPPKAkLlKJcqmR3wf03NKLZtYb2Nfd9wMuBO4sUGySg0mTYNCg0BNv1izqaEQkClkTubu/BHxaQ5M+wIhU2ylAczNrWZjwpCbvvx9WMhw5Etq1izoaEYlKIUattAIWpz3+AGhdgONKDdauhZNOgiuvhO9/P+poRCRKhZrZmVmIr/aq5uDBgzfdTyQSJBKJAn185bn00tAL//nPo45ERAopmUySTCbzek9Oo1bMrC3wpLsfXM1rdwFJd3849Xg+cLS7L81op1ErBXLffWF3n2nToGnTqKMRkWKqrwlBTwDnpD6wC7AyM4lL4cyaBb/8JTz2mJK4iARZSytm9hBwNNDCzBYD1wCNAdz9bncfZ2a9zewd4HPg/GIGXMlWrgxri99+Oxx0UNTRiEip0ISgmHAPFzfbtAmJXEQqgzaWKCN/+hMsWQKPPhp1JCJSapTIYyCZhJtvDhc3tUGEiGTS6oclbsmSsPP9yJGhrCIikkmJvIRt2ABnnAE/+lFY2VBEpDpK5CXs978PP3/722jjEJHSphp5iZowIaxqOGOGVjQUkZqpR16CPv4YzjkHRo2C3XaLOhoRKXVK5CVm/fpQF7/oorD3pohINkrkJeb3vw+llEGDoo5EROJCNfIS8txzMGyY6uIikh8l8hLx0Udw7rnw4IPQUttyiEgeVFopAevXQ79+8JOfQLduUUcjInGjRbNKwMCBMHMmjBsHDfTVKiJptGhWDDz9dBhmOGOGkriI1I4SeYTeew8GDIDRo2GXXaKORkTiSn3AiHz1FZx2GvziF9C1a9TRiEicqUYekcsvh0WLYOxYsBqrXyJSyVQjL1GPPRYS+IwZSuIiUnfqkdezBQvgyCPDCJXDDos6GhEpdbn0yFUjr0erV8PJJ8PvfqckLiKFox55PXGH/v1hm21g+HCVVEQkN6qRl5AhQ0JZZdIkJXERKayspRUz62lm881sgZn9qprXE2b2mZnNTN20bl+GiRPhxhvDePFttok6GhEpNzX2yM2sIfAXoDvwITDNzJ5w9zczmr7g7n2KFGOsLV4c1hcfNQr22ivqaESkHGXrkXcG3nH3d919HfAwcGI17VQsqMbatXDKKXDFFdC9e9TRiEi5ypbIWwGL0x5/kHounQNHmtlsMxtnZgcVMsA4u/RS2HPPMHtTRKRYsl3szGWYyQygjbuvMbNewBigXZ0ji7k77oDJk2HKFF3cFJHiypbIPwTapD1uQ+iVb+Luq9LuP2Nmd5jZTu6+IvNggwcP3nQ/kUiQSCRqEXLpGz8+bNk2eTJsv33U0YhInCSTSZLJZF7vqXEcuZk1At4CjgU+AqYC/dMvdppZS2CZu7uZdQYedfe21RyrIsaRz50bNod4/HE46qiooxGRuKvzOHJ3X29mlwDPAg2BYe7+ppldlHr9bqAv8GMzWw+sAfoVJPoYWrIETjgBbrlFSVxE6o9mdhbImjWhJ967N1xzTdTRiEi5yKVHrkReAFVVcPrpsNVWYby4Lm6KSKFoin49GTQIPv4YJkxQEheR+qdEXkdDh8Kjj8Krr0KTJlFHIyKVSIm8Dh56CAYPDmuptGgRdTQiUqmUyGtpzJgw9X7CBGhX8dOfRCRKSuS18OyzcOGF8Mwz8O1vRx2NiFQ6JfI8vfginHVW6JF36hR1NCIi2uotL1OnQt++8PDD0LVr1NGIiARK5DmaPRt+8IOwTduxx0YdjYjI15TIc/DSS9CjB/zlL2EKvohIKVEiz2L06LA5xKhRcOqpUUcjIrI5XeyswR13wLXXhmVpO3aMOhoRkeopkVfDPUy7//vfw673e+8ddUQiIlumRJ5h3Tq46KKwrvjkybDLLlFHJCJSMyXyNMuWwdlnQ6NGYdr9dttFHZGISHa62Jny/PPQoUOohY8ZoyQuIvFR8T3ydevCRhAjRoRb9+5RRyQikp+KTuTvvgv9+8OOO8LMmbDrrlFHJCKSv4osrbjDI49A585hyv1TTymJi0h8VVyPfMYMuPLKcGFz3Dj4zneijkhEpG4qpkf+4Ydw3nlhc+TTTw9rpyiJi0g5KPtEvnp1uJjZvj3svju8/TZcfHEYYigiUg7KNp299x4MGwb33guJRCip7LVX1FGJiBRe1h65mfU0s/lmtsDMfrWFNrelXp9tZh0KH2Zu1q2Dxx8P5ZOOHWHlSnjuOXjwQSVxESlfNSZyM2sI/AXoCRwE9DezAzPa9Ab2dff9gAuBO4sUa7XWrYNXXoGrrw7J+uabw5DCDz6A224r3FZsyWSyMAeKQJxjB8UfNcVf+rL1yDsD77j7u+6+DngYODGjTR9gBIC7TwGam1nLgkeasn49TJkC118PPXvCzjvDT34SEvqECWHt8LPPhm22KeznxvmXIc6xg+KPmuIvfdlq5K2AxWmPPwAOz6FNa2BpbYP68sswPHDhQli0KNw23p8/H9q2DXXviy8OZZOddqrtJ4mIxF+2RO45HsdyeV+vXlBV9c3b2rWwalW4rV4dfgK0aAH77BOWkN1nHzj++HB///1DL1xERAJz33KuNrMuwGB375l6PBCocvcb0trcBSTd/eHU4/nA0e6+NONYuX4piIhIGnfP7Cx/Q7Ye+XRgPzNrC3wEnA70z2jzBHAJ8HAq8a/MTOK5BCIiIrVTYyJ39/VmdgnwLNAQGObub5rZRanX73b3cWbW28zeAT4Hzi961CIiskmNpRURESl9RZ+in8uEolJlZsPNbKmZzYk6ltowszZmNtHM3jCzuWZ2WdQx5cPMmpjZFDObZWbzzOyPUceULzNraGYzzezJqGOpDTN718xeT/0bpkYdTz7MrLmZ/cPM3kz9/nSJOqZcmdn+qXO+8fZZTf9/i9ojT00oegvoDnwITAP6u/ubRfvQAjKz7wKrgZHufnDU8eTLzHYDdnP3WWbWFHgN+GFczj+AmW3r7mvMrBEwCbjK3SdFHVeuzOznQCdge3fvE3U8+TKzfwGd3H1F1LHky8xGAC+4+/DU78927v5Z1HHly8waEPJWE/2+AAACU0lEQVRnZ3dfXF2bYvfIc5lQVLLc/SXg06jjqC13X+Lus1L3VwNvAntEG1V+3H1N6u5WhOs0sUkoZtYa6A3cy+ZDdOMkdrGbWTPgu+4+HML1vjgm8ZTuwMItJXEofiKvbrJQqyJ/plQjNfKoAzAl2kjyY2YNzGwWYYLZRHefF3VMeRgC/AKoijqQOnBggplNN7P/jjqYPHwL+MTM7jOzGWY21My2jTqoWuoHPFhTg2Incl1JLQGpsso/gJ+leuax4e5V7n4oYbbw98wsEXFIOTGzE4Bl7j6TGPZo03R19w5AL+CnqXJjHDQCOgJ3uHtHwoi6X0cbUv7MbCvgB8Dfa2pX7ET+IdAm7XEbQq9c6omZNQYeA0a5+5io46mt1J/FTwNx2Q7kSKBPqsb8EHCMmY2MOKa8ufvHqZ+fAI8TyqVx8AHwgbtPSz3+ByGxx00v4LXU+d+iYifyTROKUt8spxMmEEk9MDMDhgHz3P2WqOPJl5m1MLPmqfvbAMcBM6ONKjfufrW7t3H3bxH+NP4/dz8n6rjyYWbbmtn2qfvbAT2AWIzgcvclwGIza5d6qjvwRoQh1VZ/QkegRkXdWGJLE4qK+ZmFZGYPAUcDO5vZYuC37n5fxGHloytwFvC6mW1MgAPdfXyEMeVjd2BE6qp9A+B+d38+4phqK45lxpbA46E/QCPgAXd/LtqQ8nIp8ECqE7mQmE1WTH15dgeyXpvQhCARkZgr+z07RUTKnRK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjM/T8bTzxBOTN6TgAAAABJRU5ErkJggg==) ## 数值积分 数值积分: $$F(x) = \lim_{n \rightarrow \infty} \sum_{i=0}^{n-1} f(x_i)(x_{i+1}-x_i) \Rightarrow F(x) = \int_{x_0}^{x_n} f(x) dx$$ 导入贝塞尔函数: In [17]: ```py from scipy.special import jv ``` In [18]: ```py def f(x): return jv(2.5, x) ``` In [19]: ```py x = np.linspace(0, 10) p = plt.plot(x, f(x), 'k-') ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczXX///HHy1K2EqlMsu+T0W9sKWFs0UxX8lOWyKipK1kGaZkuhK6ISiHL2AalqMs6XBVizkiLRk0MMxiFLFGSNNbB+/vHDM3FDDPnzJz3WV73283t9vmc8zmf99O51cv7vD+fz/stxhiUUkr5h0K2AyillHIfLfpKKeVHtOgrpZQf0aKvlFJ+RIu+Ukr5ES36SinlR1wu+iLSQUS2i0iqiLyUzfshIvKniCRm/hnmaptKKaWcU8SVD4tIYWAy0BY4ACSISKwxJuWyQ+ONMQ+50pZSSinXudrTbwLsMsbsMcakAwuBjtkcJy62o5RSKh+4WvQrAPuy7O/PfC0rA9wrIptF5BMRCXSxTaWUUk5yaXiHjIJ+Ld8DFY0xJ0XkAWAZUMvFdpVSSjnB1aJ/AKiYZb8iGb39S4wxf2XZ/lREpopIWWPM0azHiYhOAqSUUk4wxuR6CN3V4Z1NQE0RqSIi1wFdgdisB4jIbSIimdtNALm84F9kjNE/xjBixAjrGTzlj34X+l3od3H1P3nlUk/fGHNORPoDq4DCwGxjTIqIPJP5/nTgEeBZETkHnAS6udKmUkop57k6vIMx5lPg08tem55lewowxdV2lFJKuU6fyPVAISEhtiN4DP0u/qbfxd/0u3CeODMmVBBExHhKFqWU8hYignHjhVyllFJeRIu+Ukr5ES36SinlR1y+e0ep3Dhx4gQHDhzg4MGDVKlShSpVqtiOpJRf0qKv8tWZM2f48MMP2bBhAwcOHGD//v0cOHCAU6dOcccddxAQEEBqaiolSpSgdevWtG7dmlatWhEQEGA7ulJ+Qe/eUfni6NGjREdH8+6773LXXXfRqVMnKlasSIUKFbjjjjsoW7YsmQ9mY4whJSWFdevWsW7dOhwOBwEBAbRp04bnnntOfwUolQd5vXtHi75yye7du5kwYQLvv/8+Dz30EEOGDCEoKChP5zh//jybN29myZIlREdH869//YvIyEiKFNEfokpdixZ95Ra7du1i6NChrF27lqeeeooBAwZQocLls2rnXWpqKn369OHYsWPMmDGDhg0b5kNapXyX3qevCtySJUu49957CQ4OZvfu3YwdOzZfCj5AzZo1+fzzz4mMjCQ0NJQhQ4aQlpaWL+dWSmnRV3mQnp7O888/z3PPPcd///tfoqKiuOGGG/K9HREhPDycrVu38ttvv1GvXj0++eSTfG9HKX+kwzsqVw4ePEjXrl254YYbeP/997n55pvd1vaaNWt48sknGTZsGM8884zb2lXKG+jwjsp369ato1GjRnTo0IGVK1e6teADtGvXDofDwejRo4mOjnZr20r5Gr09QuXowoULjB07lnfffZf58+fTpk0ba1mqV69OXFwcrVu3BqBPnz7WsijlzbToq2wZY4iMjCQhIYFNmzbl24VaV1SvXp1169bRunVrjDE8++yztiMp5XVcHt4RkQ4isl1EUkXkpasc11hEzonI/3e1TVWwjDG8+OKLbNy4kdWrV3tEwb/oYuEfO3Ys06ZNsx1HKa/jUk9fRAoDk4G2ZCySniAiscaYlGyOGwd8BuT6goOyY9SoUaxevZq4uDhKly5tO84VLg71tGrVCmMMffv2tR1JKa/h6vBOE2CXMWYPgIgsBDoCKZcdNwBYBDR2sT1VwMaNG8dHH31EfHw8ZcuWtR0nR9WqVbtU+EVEh3qUyiVXi34FYF+W/f3A3VkPEJEKZPxD0JqMoq/3ZXqoSZMmMWPGDNavX8+tt95qO841XSz89957L/Xq1aN58+a2Iynl8Vwt+rkp4BOAKGOMkYwZt3Ic3hk5cuSl7ZCQEF0H041mzpzJ+PHjiY+P96gx/GupVq0aMTExPPbYYyQmJlKuXDnbkZQqUA6HA4fD4fTnXXo4S0SaAiONMR0y918GLhhjxmU55if+LvTlgJPA08aY2MvOpQ9nWTJ//nyioqKIi4ujZs2atuM45cUXXyQ5OZnY2FgKFdLHT5T/cOuEayJSBNgBtAEOAt8C3S+/kJvl+DnACmPMkmze06JvQVxcHN27d2fdunUEBgbajuO09PR0WrRoQefOnXn++edtx1HKbfJa9F0a3jHGnBOR/sAqoDAw2xiTIiLPZL4/3ZXzq4L1888/89hjjzF//nyvLvgARYsWZeHChTRp0oT77ruPpk2b2o6klEfSuXf81OnTp2nevDldunThhRdesB0n3yxbtoxBgwaRmJhImTJlbMdRqsDpfPrqmowxREREcOLECRYuXHhpRStfMXDgQH7++WeWLFnic383pS6nE66pa4qOjiYhIYGYmBifLIpvvPEG+/btY/LkybajKOVxtKfvZ7766is6derEl19+SY0aNWzHKTA//vgjTZs25bPPPtPVt5RP056+ytEvv/xCly5dmDNnjk8XfMiYqmHSpEn06tWL9PR023GU8hha9P3E2bNneeSRR3jmmWcIDQ21HcctunXrRsWKFZk0aZLtKEp5DB3e8RP9+/dn3759LF261K8eXkpNTeWee+5h8+bNXvWksVK5pcM76gqxsbF88sknvPfee35V8CFjofU+ffowZMgQ21GU8gja0/dxhw4dIjg4mEWLFtGsWTPbcaw4efIkd955J7NmzbK6+pdSBUF7+uqSi/fjR0RE+G3BByhRogQTJkygf//+nD171nYcpazSou/DoqOj+fXXXxkxYoTtKNY99NBDVKtWjXfeecd2FKWs0uEdH7Vjxw7uu+8+NmzYQO3atW3H8Qg//vgjd999N4mJiVSsWNF2HKXyhQ7vKNLT0+nRowevvvqqFvwsqlevTv/+/Rk8eLDtKEpZoz19HzRs2DASExNZuXKlT06z4IpTp05Rr149pk6dSvv27W3HUcplOuGan9uwYQOPPvooP/zwA7fddpvtOB7pv//9L4MHDyYpKYnrr7/edhylXKLDO37s+PHj9OrVi+nTp2vBv4qwsDDq1q3LhAkTbEdRyu20p+9DIiIiKFy4MDNmzLAdxeNt376d5s2bk5qayk033WQ7jlJOc3tPX0Q6iMh2EUkVkZeyeb+jiGwWkUQR+U5EWrvaprrSqlWrWLt2LePHj7cdxSvUqVOHsLAw3n77bdtRlHIrV9fILUzGGrltgQNAApetkSsiJY0xJzK3g4ClxpgrpnjUnr7zjh8/TlBQEDNnzuT++++3Hcdr7N69m0aNGrFjxw7KlStnO45STnF3T78JsMsYs8cYkw4sBDpmPeBiwc9UCjjiYpvqMlFRUbRt21YLfh5VrVqVrl27Mm7cONtRlHIblxZGByoA+7Ls7wfuvvwgEXkYeB0IALQy5aP4+HhiY2PZunWr7SheaejQoQQFBTF48GBuv/1223GUKnCuFv1cjccYY5YBy0SkOfA+kO0TQyNHjry0HRISQkhIiIvxfNvJkyeJiIhg6tSpejHSSRUqVOCJJ55gzJgxuryi8goOhwOHw+H0510d028KjDTGdMjcfxm4YIzJ8feyiPwINDHG/H7Z6zqmn0fPP/88Bw4cYMGCBbajeLXffvuNOnXq8N1331GlShXbcZTKE7c+nCUiRci4kNsGOAh8y5UXcqsDPxljjIg0AP5jjKmezbm06OfBxo0b6dixI0lJSdxyyy2243i9YcOGcfDgQWJiYmxHUSpP8lr0XRreMcacE5H+wCqgMDDbGJMiIs9kvj8d6Az0EpF0IA3o5kqbCs6cOcOTTz7JxIkTteDnk+eff56aNWuyc+dOatWqZTuOUgVGH87yQsOHDycpKYmlS5fq3Dr5aMyYMSQlJelwmfIqOveOj9u8eTPt2rXjhx9+0LtN8llaWho1atRg9erV1K9f33YcpXJF597xYefOnSMiIoKxY8dqwS8ApUqV4qWXXuKVV16xHUWpAqNF34tMmjSJ0qVL88QTT9iO4rOeffZZNm3aREJCgu0oShUIHd7xErt376Zx48Z888031KhxxSwWKh9NmjQJh8PBkiVLbEdR6pp0TN8HGWN44IEHCAkJISoqynYcn3fixAmqVq3KF198oSuPKY+nY/o+6MMPP+TQoUMMGTLEdhS/ULJkSfr168ebb75pO4pS+U57+h7uyJEj1KtXjxUrVtC4cWPbcfzG77//Ts2aNUlKSqJChQq24yiVIx3e8THh4eGULVuWd955x3YUvzNo0CCKFi2qPX7l0bTo+5A1a9bw1FNPsW3bNkqVKmU7jt/5+eefCQ4OZteuXZQpU8Z2HKWypWP6PuLkyZP06dOHadOmacG3pFKlSoSFhTFt2jTbUZTKN9rT91Avvvgi+/bt0ykBLNu6dStt27Zl9+7dFC9e3HYcpa6gwzs+IDExkfbt25OUlMRtt91mO47f+8c//kFoaCjPPvus7ShKXUGLvpc7d+4cTZo0ITIykt69e9uOo4ANGzYQHh7Ojh07KFLE1XWHlMpfOqbv5caPH0+5cuUIDw+3HUVluu+++wgICGDx4sW2oyjlMu3pe5DU1FTuueceEhISqFq1qu04KosVK1bwyiuv8P333+t01sqjaE/fS124cIGnn36aoUOHasH3QGFhYZw9e5Y1a9bYjqKUS1wu+iLSQUS2i0iqiLyUzfs9RGSziGwRkS9FRCcqz8asWbM4deoUkZGRtqOobBQqVIiXXnqJceNyXP5ZKa/g6hq5hclYI7ctcABI4Mo1cu8Bko0xf4pIBzIWUm+azbn8dnjn4MGD3HXXXaxbt46goCDbcVQO0tPTqV69OkuXLqVhw4a24ygFuH94pwmwyxizxxiTDiwEOmY9wBjztTHmz8zdjcAdLrbpU4wx9OvXj2effVYLvocrWrQo/fv3Z+LEibajKOU0V+8/qwDsy7K/H7j7KsdHAJ+42KZPWbx4MTt27GDhwoW2o6hceOqpp6hevTq//PILAQEBtuMolWeuFv1cj8eISCvgSaBZTseMHDny0nZISAghISEuRPN8R48eJTIykkWLFnH99dfbjqNyoWzZsnTr1o3o6GhGjRplO47yQw6HA4fDwYULF1i5cmWeP+/qmH5TMsboO2TuvwxcMMaMu+y4+sASoIMxZlcO5/K7Mf0nn3ySkiVL8u6779qOovIgJSWFkJAQ9u7dS7FixWzHUX5q8ODBpKSksGrVKreO6W8CaopIFRG5DugKxGY9QEQqkVHwe+ZU8P3R6tWrWbt2LWPGjLEdReVR3bp1CQ4O1iE5Zc3777/PihUrnJqby+WHs0TkAWACUBiYbYx5XUSeATDGTBeRWUAn4OfMj6QbY5pkcx6/6ekfO3aMoKAg5syZQ9u2bW3HUU747LPPiIqKIjExUR/WUm713Xff0aFDB+Li4qhXr57OveMNevfuTcmSJZkyZYrtKMpJFy5cIDAwkOnTp9OyZUvbcZSf+PXXX2ncuDFvv/02nTt3BvJ+y6bOHuVmsbGxbNiwgR9++MF2FOWCQoUKMXDgQCZMmKBFX7lFeno6jz76KI8//vilgu8M7em70e+//05QUBAfffQRzZs3tx1HuejEiRNUrlyZb7/9lmrVqtmOo3xcZGQkP/74I7GxsRQuXPjS6zr3jgfr168f3bp104LvI0qWLMmTTz7J5MmTbUdRPm7u3Ll89tlnfPDBB/9T8J2hPX03+fjjj3nllVdITEzUFZh8yMV1dPfs2cMNN9xgO47yQQkJCYSGhhIfH09gYOAV72tP3wMdPnyYyMhI5s2bpwXfx1SqVInWrVszd+5c21GUD0pLS6Nr165Mnz4924LvDO3pFzBjDJ06dSIwMFDvyfdRX375Jb1792bHjh0UKqT9KJV/+vXrx8mTJ5kzZ06Ox+jdOx5m/vz5/PTTT3z00Ue2o6gCcu+991K6dGk++eQTHnzwQdtxlI+Ii4tj+fLlJCUl5et5tVtSgPbu3cuQIUOYN2+ezq3jw0SEQYMGMWHCBNtRlI9IS0sjIiKC6dOnU6ZMmXw9tw7vFJD09HRatGhB586def75523HUQXs7NmzVK5cmbVr1+bb2KvyXwMGDOCvv/7K1bUifSLXQ0RFRbFlyxZWrlyp47x+4pVXXuHo0aN6C6dyicPhoGfPniQlJeWql69F3wOsWrWKiIgIEhMTueWWW2zHUW5y4MABgoKC2LNnDzfeeKPtOMoLnThxgvr16zNhwgT+8Y9/5OozesumZQcPHqR379588MEHWvD9TIUKFWjbti3vvfee7SjKS7388ss0a9Ys1wXfGdrTz0fnz5+nXbt2tGzZkhEjRtiOoyyIj4+nT58+JCcn6+ybKk/i4+N57LHH2Lp1a54u3mpP36IxY8ZgjGHYsGG2oyhLWrRoQZEiRVi3bp3tKMqLnDhxgoiICKKjo/P9bp3LaU8/n6xfv56uXbvy3Xffcfvtt9uOoyyKjo5m1apVLF261HYU5SWGDBnC4cOHmT9/fp4/6/YLuSLSgb8XUZmVzVKJdYA5QDAw1BgzPofzeG3RP3LkCMHBwcycOZMOHTrYjqMsS0tLo3LlyiQmJlKpUiXbcZSHS05OpmXLlmzbto1bb701z5936/COiBQGJgMdgECgu4jUveyw34EBwFuutOWpzp8/T3h4ON27d9eCrwAoVaoUjz/+ONHR0bajKA9njGHgwIEMGzbMqYLvDFfH9JsAu4wxe4wx6cBCoGPWA4wxvxljNgHpLrblkaKiojh16hSjR4+2HUV5kL59+zJr1ixOnz5tO4ryYMuXL+fgwYP07dvXbW26WvQrAPuy7O/PfM0vzJo1i+XLl7No0SKKFi1qO47yILVq1SI4OJj//Oc/tqMoD3X69Gmee+45Jk6c6Nb64WrR985B+HwQFxfH0KFDWblyJWXLlrUdR3mgfv366dO5Kkfjx48nODiYtm3burVdV2fZPABUzLJfkYzevlNGjhx5aTskJISQkBBnT1Wgdu7cSbdu3ViwYAG1atWyHUd5qLCwMCIjI0lISKBx48a24ygPsm/fPt5++202bdqU5886HA4cDofTbbt0946IFAF2AG2Ag8C3QHdjTEo2x44E/vL2u3eOHj1K06ZNefHFF3nqqadsx1Ee7o033iA5OVkXWVH/o3v37tSsWZNXX33V5XPZuGXzAf6+ZXO2MeZ1EXkGwBgzXUTKAwnAjcAF4C8g0BiTdtl5PL7onz17lg4dOtCgQQPeessnb0ZS+ezIkSPUqFGD1NRUnZZDARnP9PTs2ZPt27dTokQJl8+nE64VEGMM//znPzl06BDLli1zeXFi5T+eeOIJateuTVRUlO0oyrJz587RsGFDhg4dSpcuXfLlnDoNQwF56623+Pbbb/nwww+14Ks86devH9OmTeP8+fO2oyjLZs6cSZkyZXj00UetZdCinwvvvPMO06ZNY+XKldxwww224ygv06hRIwICAli5cqXtKMqi33//nREjRjBp0iSrk/Fp0b+Gt956iylTpuBwOKhYseK1P6BUNvr168eUKVNsx1AW/fvf/+aRRx6hfv36VnPomP5VvPHGG8ycOZO4uDjuuOMO23GUFzt9+jSVKlViw4YNepuvH9q9ezeNGjUiOTmZ2267LV/PrWP6+eT1119n1qxZOBwOLfjKZcWKFSMiIoKpU6fajqIsGD58OJGRkfle8J2hPf1sjB49mvfee4+4uDidJlnlm71799KgQQP27t1LqVKlbMdRbpKYmEhoaCg7d+4skGuC2tN30auvvsr8+fNxOBxa8FW+qly5Ms2bN+eDDz6wHUW5UVRUFMOHD/eYm0C06Gc6d+4cL7zwAgsXLiQuLo6AgADbkZQPunhB11N+1aqC9fnnn/PTTz/x9NNP245yiRZ94PDhw9x///1s3ryZ9evXU758eduRlI9q06YNZ86cYcOGDbajqAJ24cIFXnrpJcaMGeNRs/D6fdH/6quvaNSoEc2aNePTTz+lXLlytiMpH1aoUCG9fdNPfPTRRxQuXJhHHnnEdpT/4bcXco0xTJ48mddee42YmBjCwsLc1rbyb3/++SdVqlQhOTlZhxF91NmzZ6lTpw4xMTEFPluwXsjNhbS0NHr06EFMTAxff/21FnzlVqVLl6Zr167MnDnTdhRVQKZPn06dOnU8cnp4v+vpJyYm0rNnT+6++26mTJlC8eLFC7xNpS6XlJREhw4d2LNnj0eN9yrXHT9+nFq1arF69Wq3PH2rPf0c7N+/n969exMaGsqLL75ITEyMFnxlTVBQEDVq1GDZsmW2o6h89tZbb9G+fXvr0y3kxOeL/l9//cXw4cO56667qFChAjt27CA8PNx2LKX0gq4POnToEFOmTMmXxVEKis8W/XPnzjFjxgxq167N3r17SUxMZPTo0dx44422oykFQKdOndi5cydbt261HUXlk9GjRxMeHk7lypVtR8lRfqyc1YG/V86aZYwZl80xk4AHgJNAb2NMYjbH5MuY/m+//cZ//vMfpk6dSrly5Rg/fjwNGzZ0+bxKFYRRo0Zx6NAhpk2bZjuKctHFaTZSUlK49dZb3dauW1fOEpHCZKyR25aMRdITuGyNXBEJBfobY0JF5G5gojGmaTbncrroHz9+nKVLl7JgwQK++eYbQkNDCQ8P5/7777c6b7VS1/LLL79w55138tNPP3HTTTfZjqNcEBERQUBAAK+99ppb281r0S/iYntNgF3GmD2ZjS8EOgJZF0Z/CJgHYIzZKCI3ichtxpjDzjZ64cIF9uzZw6ZNm/j4449Zs2YNISEh9O7dm8WLF1OyZEnn/0ZKuVFAQAChoaHMnj2bIUOG2I6jnLRz505iY2NJTU21HeWaXC36FYB9Wfb3A3fn4pg7gGyLvjGG06dPc+rUKU6ePMmJEydITU1l27ZtbNu2jeTkZFJSUihbtiz169enU6dOl5YgU8obDRw4kC5dujBo0CBditNLjRgxgsGDB3vFrzVXi35ux2Mu/+mR7eeKFy/OmTNnuP766ylevDjFixenRIkSVK9encDAQFq2bEnfvn0JDAzUC7LKZzRu3Jjy5cuzcuVKOnbsaDuOyqMtW7YQFxfnNQ/buVr0DwBZ1xCsSEZP/mrH3JH52hWee+45ihQpgogQEhLikU+zKVUQIiMjmThxohZ9LzR8+HCioqLctkaCw+HA4XA4/XlXL+QWIeNCbhvgIPAtV7+Q2xSYkN8XcpXydmfPnqVq1ap89tlnBAUF2Y6jcmnjxo088sgjpKamUqxYMSsZ3PpErjHmHNAfWAUkAx8ZY1JE5BkReSbzmE+An0RkFzAd6OtKm0r5ouuuu45nn32Wd99913YUlQfDhg1j+PDh1gq+M/xu7h2lPNWvv/5K7dq12bVrFzfffLPtOOoaHA4HTz31FCkpKVbnT9K5d5TyUrfeeisdO3Zk1qxZtqOoazDGMHToUEaOHOl1E+Zp0VfKg0RGRjJlyhTOnTtnO4q6ik8//ZRjx47RvXt321HyTIu+Uh6kQYMGVK5cmeXLl9uOonJw4cIFhg0bxr///W+vfK5Ci75SHubi7ZvKMy1ZsoRChQrRqVMn21GcohdylfIw6enpVKtWjdjYWIKDg23HUVmcP3+eoKAg3n77bTp06GA7DqAXcpXyekWLFqVv3756+6YH+uCDD7j55ptp37697ShO056+Uh7oyJEj1KxZk507d3LLLbfYjqP4e7HzuXPn0qJFC9txLtGevlI+oFy5cnTu3Jno6GjbUVSmmJgYatas6VEF3xna01fKQyUnJ9O6dWv27NnjVU98+qJTp05Rs2ZNli5dSuPGjW3H+R/a01fKRwQGBtKoUSPee+8921H83rRp02jcuLHHFXxnaE9fKQ8WHx/P008/TUpKilfeE+4L/vrrL2rWrMnnn39OvXr1bMe5gvb0lfIhLVq0oEyZMsTGxtqO4rcmTpxImzZtPLLgO0N7+kp5uEWLFvH222/z1Vdf2Y7id/744w9q1arF119/TY0aNWzHyZb29JXyMZ06deLXX3/lyy+/tB3F77z55ps8/PDDHlvwnaE9faW8wNSpU1m9ejXLli2zHcVvHD58mMDAQBITE6lUqZLtODnKa09fi75SXuDkyZNUrVqV+Ph46tSpYzuOXxg0aBDGGI+fB8ltRV9EygIfAZWBPUAXY8yxbI6LAcKAX40xOa4Dp0VfqasbNWoU+/fv95oFuL3Z3r17adCgAdu2baN8+fK241yVO4v+G8ARY8wbIvISUMYYE5XNcc2BNOA9LfpKOe/IkSPUqlWL5ORkjy9E3i48PJzKlSvz6quv2o5yTe4s+tuBlsaYwyJSHnAYY7L93SkiVYAVWvSVck2/fv246aabGD16tO0oPmvLli20a9eO1NRUbrzxRttxrsmdRf8PY0yZzG0Bjl7cz+bYKmjRV8plP/74I02bNmX37t2UKlXKdhyfFBYWRvv27YmMjLQdJVfyWvSLXONka4DsfkcOzbpjjDEi4nLFHjly5KXtkJAQQkJCXD2lUj6levXqtGrVitmzZzNw4EDbcXyOw+EgJSWFJUuW2I6SI4fDgcPhcPrzrg7vhBhjDolIABCnwztKFbyEhAQeffRRUlNTvW5Rbk9mjKFp06YMHDiQxx57zHacXHPnw1mxQHjmdjigNxAr5QaNGzematWqLFy40HYUn7J48WLOnj1Lt27dbEcpUK7esvkxUIkst2yKyO3ATGNMWOZxC4CWwM3Ar8Arxpg52ZxPe/pK5VJcXBz//Oc/SUlJoUiRq47SqlxIT0/nzjvvZPLkydx///224+SJPpyllJ8ICQnhiSeeIDw8/NoHq6uKjo5m0aJFrFmzhoz7UryHFn2l/ER8fDwRERGkpKTo2L4L0tLSqFWrFrGxsTRq1Mh2nDzTCdeU8hMtW7akcuXKvP/++7ajeLUJEybQokULryz4ztCevlJebMOGDfTq1YsdO3Zob98Jv/32G3Xq1GHjxo1eO5Om9vSV8iP33XcfNWrUYO7cubajeKXXXnuN7t27e23Bd4b29JXycl9//TXdu3dn587fphsOAAANG0lEQVSdXHfddbbjeI3k5GRatGhBcnIyt956q+04TtOevlJ+5p577qFu3brExMTYjuI1jDEMGjSIYcOGeXXBd4b29JXyAd9++y2PPPIIqampXH/99bbjeLzly5fz8ssvs3nzZq+/FqI9faX8UJMmTQgKCmLWrFm2o3i806dPM3jwYCZOnOj1Bd8Z2tNXykds2rSJhx9+mF27dlGsWDHbcTzWmDFjSEhIYOnSpbaj5At9OEspP/bQQw/Rrl07BgwYYDuKR9q/fz933XUXCQkJVKtWzXacfKFFXyk/lpiYSFhYGLt27aJEiRK243icHj16ULVqVV577TXbUfKNFn2l/FzXrl0JDAxkxIgRtqN4lA0bNtC9e3e2b99OyZIlbcfJN1r0lfJzFxf1/v7776lcubLtOB7h/PnzNG7cmBdeeIHu3bvbjpOv9O4dpfxc5cqViYyM5IUXXrAdxWPMnj2bkiVL+vxc+bmhPX2lfNCpU6eoW7cuc+fO9ftlR//44w/q1q3Lp59+SnBwsO04+U6Hd5RSACxatIhXX32V77//3q8XWunbty8XLlwgOjradpQC4dbhHREpKyJrRGSniKwWkZuyOaaiiMSJyDYR2Soi3rHEvFJernPnzpQrV44ZM2bYjmJNXFwcK1asYOzYsbajeAyXevoi8gZwxBjzhoi8BJQxxkRddkx5oLwx5gcRKQV8BzxsjEm57Djt6SuVz5KSkmjTpg0pKSncfPPNtuO4VVpaGvXr1+fdd98lLCzMdpwC49bhHRHZDrQ0xhzOLO4OY0yda3xmGfCuMWbtZa9r0VeqAAwYMIDz588zdepU21Hcqn///qSlpfn8tNPuLvp/GGPKZG4LcPTifg7HVwHigTuNMWmXvadFX6kCcPToUerWrcvq1au56667bMdxC4fDQc+ePUlKSqJMmRxLkk/Ia9G/5tUdEVkDlM/mraFZd4wxRkRyrNqZQzuLgIGXF/yLRo4ceWk7JCTE7+86UCo/lC1bllGjRhEZGYnD4fC6hb/zKi0tjYiICKZPn+6TBd/hcOBwOJz+fH4M74QYYw6JSAAQl93wjogUBVYCnxpjJuRwLu3pK1VAzp8/T8OGDfnXv/5Fly5dbMcpUAMGDOD48ePMmzfPdhS3cPfwzhvA78aYcSISBdyUzYVcAeZlHjf4KufSoq9UAfriiy/o0aMHSUlJlC5d2nacAhEfH3/p7+iLvfzsuLvolwU+BioBe4AuxphjInI7MNMYEyYi9wHrgS3AxcZeNsZ8dtm5tOgrVcD69u3LsWPH+OCDD3xumOfEiRPUr1+fiRMn8uCDD9qO4zb6cJZSKkcnT56kcePGREVF8fjjj9uOk68iIyP5888//WZY5yIt+kqpq9qyZQtt2rThm2++oXr16rbj5It169bRq1cvvxrWuUgnXFNKXVX9+vUZNmwYjz32GOnp6bbjuGzPnj306NGDuXPn+l3Bd4b29JXyQ8YYwsLCCA4OZvTo0bbjOO3EiRM0a9aM3r17M2jQINtxrNDhHaVUrhw+fJjg4GA+/PBDr3wmxhhD165dKVGiBHPmzPG5C9O5pcM7Sqlcue2224iJiaFXr14cPXrUdpw8e/3119m7dy/R0dF+W/CdoT19pfzcoEGD2LdvH4sWLfKa4rly5UqeeeYZEhISuP32223HsUp7+kqpPBk7diy7du1i1qxZtqPkSkpKCk8++SSLFi3y+4LvDP9dWUEpBUCxYsVYsGABrVq1okKFCoSGhtqOlKNjx47RsWNHxo0bxz333GM7jlfS4R2lFADffPMNDz30EAsWLKBNmza241zh/PnzPPjgg9SqVYuJEyfajuMx9O4dpZTTvvjiCzp37szixYtp3ry57TiXnD59mp49e5KWlsaKFSsoWrSo7UgeQ8f0lVJOa968OQsWLKBz58588803tuMAcPz4cR544AEKFSrE8uXLteC7SIu+Uup/tGnThnnz5tGxY0e+//57q1kOHz5MSEgIgYGBLFiwgOuvv95qHl+gRV8pdYUHHniAGTNmEBoaSlJSkpUMP/30E82aNaNjx45MnjyZwoULW8nha/TuHaVUtjp27MiZM2do3749a9eupW7dum5r+4cffiAsLIzhw4fTp08ft7XrD7Snr5TKUZcuXXjzzTdp3rw5kyZN4vz58wXepsPh4P7772fixIla8AuA3r2jlLqmHTt28PTTT5Oens7s2bMJDAzM9zaOHj3KmDFjmDdvHh9//DGtWrXK9zZ8kdvu3hGRsiKyRkR2ishqEbkpm2OKichGEflBRJJF5HVn21NK2VO7dm0cDgfh4eG0bNmSUaNGcfbs2Xw595kzZxg/fjy1a9fmr7/+YsuWLVrwC5ArwztRwBpjTC1gbeb+/zDGnAZaGWP+H1AfaJW5fKJSyssUKlSIPn36kJiYyHfffUeDBg1cuq3zwoULfPjhh9SpU4f4+HjWr1/P9OnTCQgIyMfU6nJOD++IyHagpTHmsIiUBxzGmDpXOb4EEA+EG2OSs3lfh3eU8hLGGD7++GMGDRpEw4YNCQkJoUWLFgQHB1/zPvojR46wceNGRowYQaFChXjzzTdp2bKlm5L7Hrc9kSsifxhjymRuC3D04v5lxxUCvgeqA9OMMS/mcD4t+kp5mWPHjrFmzRrWr1/P+vXr2b17N02bNqVFixbcd999nDlzhpSUlP/5k56ezp133smAAQPo0qULhQrp/SSuyNeiLyJrgPLZvDUUmJe1yIvIUWNM2aucqzSwCogyxjiyed+MGDHi0n5ISIhXLuyglD87evQoX375JevXr2fDhg2UKFGCunXr/s+f8uXLe80Uzp7I4XDgcDgu7Y8aNcptPf3tQIgx5pCIBABxVxveyfzMcOCUMeatbN7Tnr5SSuWRO+feiQXCM7fDgWXZhCl38a4eESkOtAMSXWhTKaWUC1zp6ZcFPgYqAXuALsaYYyJyOzDTGBMmIvWBuWT841IIeN8Y82YO59OevlJK5ZFOrayUUn5Ep1ZWSimVIy36SinlR7ToK6WUH9Gir5RSfkSLvlJK+REt+kop5Ue06CullB/Roq+UUn5Ei75SSvkRLfpKKeVHtOgrpZQf0aKvlFJ+RIu+Ukr5ES36SinlR7ToK6WUH3G66ItIWRFZIyI7RWT1xRWycji2sIgkisgKZ9tTSinlOld6+lHAGmNMLWBt5n5OBgLJgK6SkgtZFz32d/pd/E2/i7/pd+E8V4r+Q8C8zO15wMPZHSQidwChwCwg16u7+DP9D/pv+l38Tb+Lv+l34TxXiv5txpjDmduHgdtyOO4d4AXgggttKaWUygdFrvamiKwBymfz1tCsO8YYIyJXDN2IyIPAr8aYRBEJcSWoUkop1zm9MLqIbAdCjDGHRCQAiDPG1LnsmDHA48A5oBhwI7DYGNMrm/PpeL9SSjkhLwuju1L03wB+N8aME5Eo4CZjTI4Xc0WkJfC8MeYfTjWolFLKZa6M6Y8F2onITqB15j4icruI/DeHz2hvXimlLHK6p6+UUsr7WH8iV0Q6iMh2EUkVkZds57FFRCqKSJyIbBORrSISaTuTbfpQXwYRuUlEFolIiogki0hT25lsEZGXM/8fSRKRD0XketuZ3EVEYkTksIgkZXkt1w/JXmS16ItIYWAy0AEIBLqLSF2bmSxKBwYbY+4EmgL9/Pi7uEgf6sswEfjEGFMXqA+kWM5jhYhUAZ4GGhhjgoDCQDebmdxsDhm1Mqu8PCQL2O/pNwF2GWP2GGPSgYVAR8uZrDDGHDLG/JC5nUbG/9i3201ljz7Ul0FESgPNjTExAMaYc8aYPy3HsuU4GZ2jEiJSBCgBHLAbyX2MMV8Af1z2cq4eks3KdtGvAOzLsr8/8zW/ltmjCQY22k1ilT7Ul6Eq8JuIzBGR70VkpoiUsB3KBmPMUWA88DNwEDhmjPncbirrcvuQ7CW2i76//2y/goiUAhYBAzN7/H4n60N9+HEvP1MRoAEw1RjTADhBLn7C+yIRqQ4MAqqQ8Su4lIj0sBrKg5iMu3KuWVNtF/0DQMUs+xXJ6O37JREpCiwG5htjltnOY9G9wEMishtYALQWkfcsZ7JlP7DfGJOQub+IjH8E/FEj4CtjzO/GmHPAEjL+W/Fnh0WkPEDmQ7K/XusDtov+JqCmiFQRkeuArkCs5UxWiIgAs4FkY8wE23lsMsb8yxhT0RhTlYwLdeuye4rbHxhjDgH7RKRW5kttgW0WI9m0HWgqIsUz/39pS8aFfn8WC4RnbocD1+wsXnXunYJmjDknIv2BVWRciZ9tjPHLOxOAZkBPYIuIJGa+9rIx5jOLmTyFvw8DDgA+yOwY/Qg8YTmPFcaYzZm/+DaRca3ne2CG3VTuIyILgJZAORHZB7xCxkOxH4tIBLAH6HLN8+jDWUop5T9sD+8opZRyIy36SinlR7ToK6WUH9Gir5RSfkSLvlJK+REt+kop5Ue06CullB/Roq+UUn7k/wDY8//eOSbGSgAAAABJRU5ErkJggg==) ### `quad` 函数 Quadrature 积分的原理参见: [http://en.wikipedia.org/wiki/Numerical_integration#Quadrature_rules_based_on_interpolating_functions](http://en.wikipedia.org/wiki/Numerical_integration#Quadrature_rules_based_on_interpolating_functions) quad 返回一个 (积分值,误差) 组成的元组: In [20]: ```py from scipy.integrate import quad interval = [0, 6.5] value, max_err = quad(f, *interval) ``` 积分值: In [21]: ```py print value ``` ```py 1.28474297234 ``` 最大误差: In [22]: ```py print max_err ``` ```py 2.34181853668e-09 ``` 积分区间图示,蓝色为正,红色为负: In [23]: ```py print "integral = {:.9f}".format(value) print "upper bound on error: {:.2e}".format(max_err) x = np.linspace(0, 10, 100) p = plt.plot(x, f(x), 'k-') x = np.linspace(0, 6.5, 45) p = plt.fill_between(x, f(x), where=f(x)>0, color="blue") p = plt.fill_between(x, f(x), where=f(x)<0, color="red", interpolate=True) ``` ```py integral = 1.284742972 upper bound on error: 2.34e-09 ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjXX/x/HXZ86ItFhyk2RpoeSmaCH93KZIIyVSku6yVUokJamkkbuivUQphe7QpiS7W52QIpFkxp617Et2M+f6/P6YqSZmzDhn5nzP8nk+Hh6d65xreXcefOYz3+u6vpeoKsYYY+JDgusAxhhjwseKvjHGxBEr+sYYE0es6BtjTByxom+MMXHEir4xxsSRkIu+iCSLyFIRWSEij+TweZKI7BaRhVl/+oR6TGOMMcFJDGVjEfEBrwONgY3A9yIyXlXTjlj1a1VtHsqxjDHGhC7UTv8yYKWqrlHVdOAD4IYc1pMQj2OMMaYAhFr0KwDrsy1vyHovOwXqi8giEZkkIheEeExjjDFBCml4h8yCnpcFQEVV3S8iTYFxQLUQj2uMMSYIoRb9jUDFbMsVyez2/6Sqe7K9niwiQ0SktKruyL6eiNgkQMYYEwRVzfcQeqjDO/OBqiJSRUROAG4BxmdfQUTKiYhkvb4MkCML/h9U1f6o8uSTTzrPECl/7Luw78K+i2P/OV4hdfqqmiEiXYGpgA94R1XTRKRz1udDgZuAe0UkA9gPtAnlmMYYY4IX6vAOqjoZmHzEe0OzvR4MDA71OMYYY0Jnd+RGoKSkJNcRIoZ9F3+x7+Iv9l0ET4IZEyoMIqKRksUYY6KFiKBhPJFrjDEmiljRN8aYOGJF3xhj4kjIV+8YkxtVJS1tG2vXbmHLlq3s2LGTE08sxT/+UZ5//as8ZcqcQtYtHMaYMLGibwrUqlW/MnjwVCZM+JJVq77E8w4gUo6EhH8gUgrVnXjeb6j+RpEiZfjnP5O54YZk7r77KsqXP9V1fGNinl29Y0KmqrzxxiyefHIQ27bNIDGxCRkZVwFXAeeQ8ySrCqQiMoWEhKkEAt9z4YW3M3z4Q9SuXTms+Y2JZsd79Y4VfROSN96YRq9ej7Bv335UuwF3AMF07L/i871KIDCMCy64ljFj/kOtWlb8jcmLFX0TFl999Qtt2z7I5s0/ofoCmY9RKIjrAnaTmPgqGRmDaNbsST77rAtFitj1Bsbkxoq+KVSBgMdNN73IuHEDEOmBak+gWCEcaSkJCXdy0kkwffq71K1rs3EbkxO7OcsUmtWrt3LGGc0YP34c8AOqfSicgg9wPp43k/37b+Hyy/+P116bUkjHMSa+WNE3+TJs2NdUrVqb7dsvwvP8QJUwHDWBQKAbqp/RvXsHWrd+JaipZI0xf7HhHZOndu1G8N57vYGRwDWOUqxBpDnnnluPJUveoEgRn6McxkQWG9M3BUZVadToP/j975I5g/b5jhPtISGhJVWqlGfZshEkJlrhN8aKvikQ6ekZXHzxvfz88wJUJwKnu46UZT8JCddRtWollix5F5/PRihNfAv7iVwRSRaRpSKyQkQeOcZ6l4pIhojcGOoxTeE6fDiD6tVvZ8mSNaj6iZyCD1Acz/uCFSvWULv2XXie5zqQMVElpKIvIj7gdSAZuAC4VUSq57LeQGAKOd+eaSLE4cMBzjuvA7/8sg3PGw+c4jpSDk7C8yawZMlyGjR42HUYY6JKqJ3+ZcBKVV2jqunAB2TepXOkbsAnwNYQj2cKUUaGR9WqnVi3biOe9zlwoutIx3Aynvc53377Bbfc8rbrMMZEjVCLfgVgfbblDVnv/UlEKpD5g+CNrLds4D4CqSq1a3dhw4Zf8LwvgOKuI+VDaVQn8NFHfXjppa9chzEmKoQ6y2Z+CvgrQG9VVcmcRzfX4Z2UlJQ/XyclJdlzMMOoWbNnWLJkLqpfAye5jnMcqgEf0LNnGy69dBYNGtiduya2+f1+/H5/0NuHdPWOiNQDUlQ1OWv5UcBT1YHZ1lnNX4W+DLAfuEtVxx+xL7t6x5GuXd9jyJC+qM4BznAdJygJCW9TpMhLbNo0n5Ilo+mHljGhCeslmyKSCCwDGgG/AvOAW1U1LZf1hwNfqOqnOXxmRd+BV16ZTo8e/wa+IvNcfPTy+dpTtaqQljbcdRRjwiasl2yqagbQFZgKpAIfqmqaiHQWkc6h7NsUvokT03jwwduAj4n2gg8QCLzOsmXf0q3b+66jGBOx7OasOLV27S7OOacugUBvoIPrOAVoEdCYqVPn0KRJVddhjCl0dkeuyVN6eoDy5a9n585z8LxBruMUgiEULTqMbdu+5eSTi7oOY0yhsqmVTZ4aNerLzp378byXXEcpJPeSkVGRa6/9j+sgxkQcK/px5oknxjJ79ig872OgiOs4hUQIBN5k1qyhjB270HUYYyKKDe/EkZkzV5GUVA/VScClruOEwUiKFXuZHTvmceKJJ7gOY0yhsOEdk6Pduw9xzTWtgSeIj4IPcAfp6WfQosUA10GMiRjW6ceJ6tW7sXz5RjxvLPE1590GoDaTJn1J06Y1XYcxpsBZp2+O8vDDn7Bs2UQ8713iq+ADnInIf2jd+h4CAZuG2Rgr+jHu++/X8eKLXVD9ECjpOo4Tqnexf38GHTu+5zqKMc7Z8E4MS08PUKZMI/buTcbzeruO49h8RK5j1ao0zjqrlOswxhQYG94xf2re/Hn27lU8zx40Apcg0oLrr+/rOogxTlmnH6NGj/6B225rCswHKrmOEyG2Axfw6adTadnyItdhjCkQNg2DYdu2/ZQvX4eMjBSgjes4EUVkGKecMpxdu2aT+XgHY6KbDe8YGjbsjerFWME/mmpH9u49SJ8+n7iOYowT1unHmBdf/JKHH74D1cWAnbDM2ZckJt7Fzp2pNiGbiXrW6cexjRt/p1evjqi+jRX8Y7kK1ercfvsbea9qTIyxTj+GVKt2F6tXC4HAW66jRIFURJJYuXIZZ59tPyBN9Ap7py8iySKyVERWiMgjOXx+g4gsEpGFIvKDiFwV6jHN0Z5+ehIrV/6PQOBF11GixAWItKRly6ddBzEmrEJ9Rq6PzGfkNgY2At9zxDNyReQkVd2X9bom8JmqnpvDvqzTD9L69bupUuWfeN5IwH6m5t8moAazZ8/niivOch3GmKCEu9O/DFipqmtUNR34ALgh+wp/FPwsJwPbQjymOcJVV/UEmmEF/3idTkLCfdxxRz/XQYwJm1CLfgVgfbblDVnv/Y2ItBCRNGAycH+IxzTZPPvsdFatmobnPec6SlTyvIdYvXoi06YtdR3FmLBIDHH7fI3HqOo4YJyINAD+C5yX03opKSl/vk5KSiIpKSnEeLHt11/30KfPXagOBU51HSdKlSAh4SE6dUph/foPXIcxJk9+vx+/3x/09qGO6dcDUlQ1OWv5UcBT1YHH2GYVcJmqbj/ifRvTP041atzHsmX7CQSGu44S5fYB5zJ27BRuvPFC12GMOS7hHtOfD1QVkSoicgJwCzD+iEDnSNb97iJSB+DIgm+O35Ahs0hNHUcgEKsPNw+nkxDpzb332mRsJvaFVPRVNQPoCkwFUoEPVTVNRDqLSOes1VoBi0VkIfAqNjdAyHbvPkj37ncCg7CbsAqGame2bl3A++/Pcx3FmEJlN2dFofr1H2fevKUEAmNdR4kxb1G69Kds3z7FdRBj8s2mYYhxH3+8iG+/fZtA4HXXUWJQe3buTOXdd63bN7HLOv0ocuhQBiVLXs7Bg/cAnVzHiVGDKV16Ktu3j897VWMigHX6MaxVq1c5fPhUoKPrKDGsEzt2/MCYMQtdBzGmUFinHyVmzvyFhg0vBb4DjprFwhSolylXbjabNtk5ExP57MlZMcjzlNNOS2b37qtQPWpOO1Pg9gNnM378/7j++n+6DmPMMdnwTgy6995R/P77ZlQfdB0lThRHpAddu9oMnCb2WKcf4ZYt20b16v9EdQJwies4cWQPcDZ+/xwaNqzqOowxubLhnRhTufIdbNhwGp73susocSchoS/Vqm0mLW2o6yjG5MqGd2LIs89OY/36mXhef9dR4pLndWPp0o9ZvPg311GMKTDW6Ueobdv2c/rp/yQQGAw0dR0nbiUkdOPii09i3rwBrqMYkyMb3okRF1/ci0WLNhAIjHYdJc6tAS5mzZrVVK5cwnUYY45iwzsxYPTohSxYMJJA4BXXUQxV8PmSufPON10HMaZAWKcfYQ4ezKBkybocOtQNaO86jgHgJ0SS2blzNSVKFHMdxpi/sU4/yt1446ukp5cE2rmOYv5Ui4SEi+je/b+ugxgTMuv0I8iXX66mUaPLgLnAOa7jmL/5kiJF7uPAgSX4fNYrmchhnX6U8jzlxhs7I/IIVvAj0ZUEAsV45pnJroMYE5KQi76IJIvIUhFZIZkV68jPbxORRSLyk4h8IyK1Qj1mLLrzzpHs2bMd1R6uo5gcCZ73EC+8YI+nNNEt1Aej+4BlQGNgI/A9cKuqpmVb53IgVVV3i0gymQ9Sr5fDvuJ2eOfnnzdTq1ZNVKcAdVzHMbk6DJzNhx9OoHXri1yHMQYI//DOZcBKVV2jqunAB8AN2VdQ1W9VdXfW4lzgzBCPGXOaNLkfkQ5YwY90JyDSjYcfftF1EGOCFmrRrwCsz7a8Ieu93HQCJoV4zJjyyCPj2LTpRzwvxXUUkw+qd7Nu3UTmz9/gOooxQUkMcft8j8eIyJVkPvLpitzWSUlJ+fN1UlISSUlJIUSLfKtX7+T55+9D9QPgRNdxTL6Uwue7nS5dBjFv3kDXYUwc8vv9+P1+AoEAEydOPO7tQx3Tr0fmGH1y1vKjgKeqA49YrxbwKZCsqitz2VfcjemfdVZH1q0rjufZQ86jy2rgMjZvXkvZsie5DmPikKpy1113sWnTJiZOnBjWMf35QFURqSIiJwC3AH97orSIVCKz4P87t4Ifj/r3n8batV/iec+6jmKO29n4fA3o2vU910FMnHrhhReYP38+Y8aMOe5tQ745S0SaAq8APuAdVX1WRDoDqOpQERkGtATWZW2SrqqX5bCfuOn0N27cQ6VKNfG8ocA1ruOYoPhJTLyXAweWkJhot7uY8Pnss8/o1q0b3377LRUrVrRZNqPBeefdw6pVGQQCw1xHMUFTEhJq06/fAPr0SXYdxsSJH3/8kauvvpopU6Zw8cUXA3ZHbsR77rnprFgxiUDALvuLboLndefFF20mVBMeu3fv5qabbmLQoEF/FvxgWKcfRuvX/06VKjXxvLeBJq7jmJAdBCozaZKfpk2ruw5jYpiq0qpVK8qXL8/gwYP/9pkN70Sws8++i7VrBc97y3UUU0BE+nLeeVtJS3vDdRQTw15++WVGjx7N7NmzKVq06N8+s6IfoZ54YgpPP90Z1cXAqa7jmALzG3ABv/yyiipVSrsOY2LQt99+S4sWLfjuu+8466yzjvrcxvQj0NKl23n66U6ojsAKfqwpj8/XjG7d3nUdxMSg33//nbZt2/LWW2/lWPCDYZ1+IfM85YwzWrN1a0U8z2ZojE1z8fnacODASooU8bkOY2JIx44dSUxM5K23ch8Stk4/wtxzzyi2bEnF855xHcUUmrqo/oOnnjr+W+KNyc24ceOYOXMmL71UsM2idfqF6Ntv11G//sXANKC26zimUL1PiRIj2bVruusgJgZs3ryZCy+8kE8//ZT69esfc13r9CPEoUMBmjS5A5EHsYIfD25m9+7FfPFFWt6rGnMMqsqdd95Jp06d8iz4wbCiX0iuvvpZ9u8XVHu5jmLCoigid9Orl02eZ0IzatQo1q1bx5NPPlko+7fhnULw+utz6NbtRuAHjv14ARNbfgX+ydq1v1CpUgnXYUwU2rp1KzVr1mTChAlccskl+drGrtN3bNWqXVSrVhvPexVo7jqOCTOfrw3Nm1/Op592dx3FRKHbbruN8uXL88ILL+R7Gyv6DnmeUr58G7ZtK4vnDXIdxzjxDT5fBw4dWorPZ6OnJv8mTZpE165dWbx4MSedlP/nNNiJXIduumkI27Ytw/Oedx3FOFMf1eI895xdxWPyb8+ePdx7770MHTr0uAp+MKzTLyAjR86lffvrgTnAua7jGKeGUabMeLZuHZ/3qsYADz74INu3b2fkyJHHvW3YO30RSRaRpSKyQkQeyeHz80XkWxE5KCIPhXq8SLRs2XY6dmwNvIUVfANt2bZtDrNm/eI6iIkCixcv5r///S/PPx+eEYJQn5HrA5YBjYGNwPfAraqalm2dfwCVgRbATlXNcSL5aO30Dx/2KFu2GXv2/NOGdcyfEhIeonZtH/PnP+c6iolgqkrDhg1p06YNXbp0CWof4e70LwNWquoaVU0HPgBuyL6Cqm5V1flAeojHikhXXJHCnj37bJoF8zeedy8//DCcnTsPuI5iItioUaPYt28fnTt3DtsxQy36FYD12ZY3EEcXpnfr9jE//DASz/sYKOI6joko5+LzXUrPnh+4DmIi1O7du+nVqxdDhgzB5wvfRH2hFv3oG48pIGPGLOT117ugOg4o5zqOiUCBwH2MGvU6nhe3/0zMMTz55JM0a9aMunXrhvW4iSFuvxGomG25IpndflBSUlL+fJ2UlERSUlKwuypUS5Zs5t//bgEMwebVMblryuHD3Rg+fC6dOtVzHcZEkNTUVEaNGkVqaupxb+v3+/H7/UEfO9QTuYlknshtROY96PM44kRutnVTgD3RfiJ327YDVKrUmIMHr0K1v+s4JuI9T6VKi1m79j3XQUyEUFWaNm1KcnIyDzzwQMj7C/sduSLSFHgF8AHvqOqzItIZQFWHisjpZF7VcyrgAXuAC1R17xH7ifiif+hQgDPPvJkdO4rhee9j97aZvG0HziE1dQXVq//DdRgTASZNmkSPHj1YvHgxJ5xwQsj7s2kYConnKdWr38/KlUvwvMlA0Ty3MQbA5+vAVVedx7RpvV1HMY6lp6dTs2ZNXnzxRZo1a1Yg+7RpGApJ48bPs2LF13jeZ1jBN8cjELiPGTPe5PDhgOsoxrEhQ4ZQpUoVrr32WmcZrOjnw7///Q5+/2BUJwE2Za45XpcA5ejXzx6nGM+2b9/O008/zUsvvYRIvhvzAmfDO3m49973efPN3sBXQFXXcUzUeo8SJUaxa9dU10GMIw888ADp6ekMHjy4QPdrY/oF6N57P+LNNx8A/gdc4DqOiWoHgUpMnfoNTZpY8xBvVq1aRd26dUlNTaVs2bIFum8r+gXknns+ZejQLmQ+1LyW6zgmBoj0pmbNwyxa9JLrKCbMWrduzYUXXsjjjz9e4Pu2ol8A2rZ9lzFjHgcmYTdfmYKzBriYzZvXUbZs4c6ZbiLH3LlzadWqFcuXL6d48eIFvn+7eidE1133Ah988BTwNVbwTcGqgs93BQ8+ONp1EBMmqkrPnj156qmnCqXgB8OKfpbDhz3q1OnNpEnvoDoLqOY6kolBgUBXPv54sM3HEyc+//xzdu/eTbt27VxH+ZMVfWD9+r2UL38TixbNyir4FfPcxpjgNCY9fT9Dh85xHcQUsoyMDB599FEGDhwY1lk08xL3RX/GjDWcfXZ9du0qhed9CZRxHcnEtARUu/D00wV72Z6JPCNHjqRcuXIkJye7jvI3cX0i9z//mUrfvu1QfRS4H3B3w4SJJ7uAs1i0aCm1atm03LHowIEDVKtWjY8//ph69Qp3hlU7kZsP27cfpFatHvTteyeqY4DuWME34VOShISbuf/+t1wHMYVk8ODBXHLJJYVe8IMRd53+8OGLufvuf+N5VfG8t4DShX5MY472EyLXsnfvLxQvbk9diyW7du2iWrVq+P1+Lrig8G/qtE4/F7/9todatR6mY8eryMjonvWIQyv4xpVaiJzNE0987jqIKWDPP/881113XVgKfjBivtNPT1e6dv2IYcN6ItKYQGAgULC3QRsTnI84+eQh7Nnjdx3EFJBNmzZRo0YNFi5cSKVKlcJyTOv0sxw+7NG166cUL16bYcOew/PGEAgMxwq+iRwt2bt3BWPHLnYdxBSQZ555hnbt2oWt4AejIJ6clcxfT84apqoDc1jnNaApsB9or6oLc1inQDr9LVv28dhjnzBy5It43gl43pPAddiJWhOJRJ7i3HN/ZfnyN11HMSFau3YtderUIS0trcAnVTuWsM69IyI+Mp+R25jMh6R/zxHPyBWRa4GuqnqtiNQFXlXVo05ph1L0Dx3yGDLke157bQRr1nyIz1efQKALmT9nrNibSLYJqM4vv/xClSolXYcxIejYsSMVKlSgf//wPjv7eIt+YojHuwxYqaprsg7+AXADkP3B6M2BkQCqOldESopIOVXdHOxBVeG77zbzwQezmThxIqtXT0KkJJ53G/ATgcCZQf8PGRNep+PzNaV79xF8/nnoD8k2bixdupQvvviCFStWuI6Sp1CLfgVgfbblDUDdfKxzJnDMon/w4GFWrdrN8uU7SU1dz7Jla1i9ei1paYvZuXM+qvvw+S4jELgWeBzVc0L8XzHGjUCgGxMn3kFGxv0kJsbsabaY1rdvX3r27EnJkpH/21qoRT+/4zFH/uqR43annXYagUCAQ4cOcehQBqolgJJkzoVTJetPG+BFTjzxLKePHDOmoHhePcoc/J0pp57CdScU4jX7lSvDokWFt/84tXDhQmbPns3w4cNdR8mXUIv+Rv4+O1lFMjv5Y61zZtZ7R+nUqRMigs/no0GDq7nwwitDjGdMNBDG92rO4DHvct2B/YV3mGXLCm/fcaxPnz489thjnHRSeJ6R4Pf78fv9QW8f6oncRDJP5DYCfgXmcewTufWAVwr6RK4x0e7g7t1ULlWKmaqcV1gHKVoUDh4srL3HpTlz5nDrrbeyfPlyihYt6iRDWK/TV9UMoCswFUgFPlTVNBHpLCKds9aZBKwWkZXAUKBLKMc0JhYVK1GCuy+/nNcTbEw/mvTp04e+ffs6K/jBiPk7co2JFhsXLKDmxRfzC1CiMA5gnX6B+vLLL+ncuTNpaWkkJoY6Uh48uyPXmChVoU4dmlSowAi7QCHiqSqPP/44/fr1c1rwg2FF35gIcn9KCoMAz3UQc0yTJk1i7969tGnTxnWU42ZF35gIcnmnTpQsWpRJroOYXHmeR58+fXjqqadIiMJzMNGX2JgYJiJ079CBVyPomarm7z799FN8Ph8tWrRwHSUodiLXmAhzeN8+qpxyClNVqVmQO7YTuSELBALUrFmTl156KWKefWsnco2JciecdBJdrrySV6zbjzijR4+mdOnSXHPNNa6jBM06fWMi0LalS6lavTrLKMAnQFinH5L09HTOP/983n33XRo2bOg6zp+s0zcmBpQ5/3xuPucc3rTLNyPGiBEjOPvssyOq4AfDOn1jItSS8eNpfMMNrAEK5H5P6/SDdvDgQapVq8bHH39M3bpHTiTslnX6xsSIGs2bU6tUKca4DmIYOnQoF110UcQV/GBY0TcmgvXo1YuXExLyPYe5KXj79u1jwIABYX8iVmGxom9MBLumVy8CiYn8z3WQOPbaa6+RlJTEhRde6DpKgbAxfWMi3IjOnRn9zjtMCwRC25GN6R+3Xbt2UbVqVb755huqVavmOk6Owvpg9IJkRd+YnB3au5ezTz2VSaqE1Gta0T9uffr04bfffuOdd95xHSVXVvSNiUEDr72Wn6dN47+hdPtW9I/Lli1bqF69OgsWLKBy5cqu4+TKir4xMWjX2rWcXaUKi/j7s0ePixX949KjRw8yMjIYNGiQ6yjHFLaiLyKlgQ+BysAaoLWq7sphvXeBZsAWVc11KhEr+sYc24MXXUTC4sW84AU58bIV/Xxbu3YtderUITU1lXLlyrmOc0zhvE6/NzBdVasBM7KWczIciIyZiYyJYg+88QbDPY+jOitT4Pr160eXLl0ivuAHI5ROfynQUFU3i8jpgF9Vz89l3SrAF9bpGxOa2ytVosaGDfQO5t+Kdfr5kpqaSlJSEitWrKBEiUJ5cGWBCmenX05VN2e93gzE3o9EYyJM79df5xVV9rsOEsP69OlDr169oqLgB+OYD3cUkenA6Tl89Hj2BVVVEQm5TU9JSfnzdVJSEklJSaHu0piYUqN5c+qVLcu7W7bQ1XWYGDRv3jy+//57Ro0a5TpKrvx+P36/P+jtQx3eSVLVTSJSHvjKhneMKXxzR46kdYcOrFDlhOPZ0IZ3jklVady4Mbfeeit33nmn6zj5Fs7hnfFAu6zX7YBxIezLGJNPddu1o+qppzLadZAYM3XqVDZu3Ej79u1dRylUoRT9AcDVIrIcuCprGRE5Q0Qm/rGSiIwB5gDVRGS9iHQIJbAxBh7r148BIoQ4MYPJEggE6NWrFwMGDCAx8Zij3lHPbs4yJgqpKpeffDIP799Pq/xuZMM7uRo5ciRvv/02s2bNQqLswTU2n74xcUBEeKxXL/onJBDkrVomy4EDB3jiiSd47rnnoq7gB8OKvjFR6vonniChSBE+dx0kyg0aNIhLL72U+vXru44SFja8Y0wUG//YY/QdOJAFnpd3B2fDO0fZvn07559/PrNnz+a8885zHScoNuGaMXFEAwEuKV6cPocP0zKvla3oH6V79+5kZGQwePBg11GCZmP6xsQR8flI6dWLfja2f9yWLVvG6NGj/3ZTaDywTt+YKKeex6UnncRjBw9y47FWtE7/b66//noaNmxIz549XUcJiXX6xsQZSUggpU8fUkSs28+n//3vf6SlpdGtWzfXUcLOOn1jYoCqUu+UU3hg3z5uzW0l6/SBzBuxateuTUpKCjfeeMzfjaKCdfrGxCERYcBzz9FHhMOuw0S4d955h1KlStGyZZ6nvmOSdfrGxJDkMmW4fscO7svp35J1+uzYsYPq1aszefJk6tSp4zpOgbBLNo2JYws/+ohrb7mFFcDJR35oRZ+uXbsSCAR44403XEcpMFb0jYlzt555JjV++40+Rz5LN86L/o8//sg111xDamoqp512mus4BcbG9I2Jc/1HjuQVz2Ob6yARRFXp2rUr/fv3j6mCHwwr+sbEmHMbNaJNjRr08/lcR4kYo0aN4uDBg3Tq1Ml1FOdseMeYGLRt1SouqFqVr1Sp8cebcTq8s2vXLmrUqMHYsWM2s6VIAAALl0lEQVSpV6+e6zgFzoZ3jDGUOecc+tx8Mz0SEoj3VurRRx/l+uuvj8mCH4yQOn0RKQ18CFQG1gCtVXXXEetUBN4DygIKvKWqr+WwL+v0jSlA6QcPcmGJEgw4fJjmEJed/uzZs7nllltYsmQJJUuWdB2nUIS70+8NTFfVasCMrOUjpQM9VLUGUA+4T0Sqh3hcY0weihQrxsvPPMNDIhxyHcaBQ4cOcffdd/Pqq6/GbMEPRqid/lKgoapuFpHTAb+qnp/HNuOAQao644j3rdM3phBcf/rpNNi6lV5FisRVp9+/f3++//57Pv/885h+IlZYr9MXkZ2qWirrtQA7/ljOZf0qwNdADVXde8RnVvSNKQQrZs3i8n/9iwUnnEClQ/HR8y9btowrrriChQsXUrFiRddxCtXxFv08H/suItOB03P46PHsC6qqIpJr1RaRk4FPgO5HFvw/ZJ/XOikpiaSkpLziGWPyULVBA7onJdFl5ky+UI3prhcgIyOD9u3bk5KSEpMF3+/34/f7g96+IIZ3klR1k4iUB77KaXhHRIoAE4DJqvpKLvuyTt+YQnL499+pfd55pLz2GjfffLPrOIXq2WefZcaMGUybNo2EhNi/QDHcwzvPAdtVdaCI9AZKqmrvI9YRYGTWej2OsS8r+sYUom+++YbWrVvH9JUsixYtonHjxvzwww9UqlTJdZywCHfRLw18BFQi2yWbInIG8LaqNhOR/wNmAj/Bn5cMP6qqU47YlxV9YwrZPffcA8Cbb77pOEnBO3ToEJdddhk9evSgffv2ruOEjU24ZozJ1R93p77//vtceeWVruMUqEcffZTU1FTGjRsX8+ctsrOib4w5pilTpnD33XezaNEiSpXK9WK7qDJ9+nTat2/PggULKFeunOs4YWVF3xiTp27durF161bGjBkT9V3xr7/+yiWXXMKoUaNi7reX/LC5d4wxeRo4cCA//fQTo0ePdh0lJBkZGbRt25Z77rknLgt+MKzTNyZOLVy4kCZNmjB//nwqV67sOk5Q+vbty5w5c5g6dSq+OJ1K2oZ3jDH59vzzz/PJJ58wc+ZMihYt6jrOcZkwYQKdO3eOy3H87KzoG2PyTVW5+eabKVmyJG+//XbUjO//9NNPNGrUiAkTJlC3bl3XcZyyMX1jTL6JCCNGjOC7776Lmmv3N23aRPPmzXn99dfjvuAHwzp9YwwrV67kiiuuYOzYsfzf//2f6zi5OnDgAFdeeSVNmzblySefdB0nItjwjjEmKFOmTKFDhw58/fXXVKtWzXWco6Snp9O6dWuKFSvG6NGjo2YoqrDZ8I4xJijJycn079+fJk2asH79etdx/iYjI4Pbb7+dw4cPM2LECCv4IchzamVjTPy488472b17N1dffTUzZ86kbNmyriPheR4dO3Zk+/btfPHFF1F3lVGksaJvjPmbhx56iF27dpGcnMyMGTOcTtUQCATo3Lkza9euZfLkyRQrVsxZllhhwzvGmKM89dRTXHnllTRo0IB169Y5ybB3715atGjBmjVrmDBhAsWLF3eSI9ZY0TfGHEVEeOGFF+jQoQP169fnxx9/DOvxf/31V/71r39RtmxZJk+ezCmnnBLW48cyK/rGmByJCA899BAvvfQSTZo0YfLkyWE57rx587j88stp1aoVw4YNo0iRImE5brywSzaNMXmaNWsWbdu2pWXLlgwYMKBQhlrS09Pp378/Q4cOZciQIbRq1arAjxGLwnbJpoiUFpHpIrJcRKaJyFHPXxORYiIyV0R+FJFUEXk22OMZY9xp0KABP/30E9u3b6d27drMnTu3QPe/ePFi6tWrxw8//MCPP/5oBb8QhTK80xuYrqrVgBlZy3+jqgeBK1X1IqAWcGXW4xONMVGmVKlSjBo1iqeffpobbriBW2+9lZ9//jmkfaalpdG2bVsaNWrEPffcw4QJEyhfvnwBJTY5CaXoNyfzgedk/bdFTiup6v6slycAPmBHCMc0xjh20003sWLFCi666CIaN25Mq1atmD59OocOHcrX9gcOHGDChAm0adOGhg0bUqtWLVatWsVdd91lN12FQdBj+iKyU1VLZb0WYMcfy0eslwAsAM4B3lDVXrnsz8b0jYky+/btY9iwYXz44YcsWbKEpKQkkpKSqFChAmXLlqVMmTLs3LmT9evXs379eubOncuMGTOoXbs2LVu2pGPHjnZlTogKdO4dEZkOnJ7DR48DI7MXeRHZoaqlj7GvEsBUoLeq+nP4XLNPoPTHXx5jTHTYtm0b06ZNY86cOWzevJktW7awdetWSpUqRcWKFalYsSI1a9akWbNmnHbaaa7jRi2/34/f7/9zuV+/fuGZcE1ElgJJqrpJRMoDX6nq+Xls8wRwQFVfyOEz6/SNMeY4hXPCtfFAu6zX7YBxOYQp88dVPSJyInA1sDCEYxpjjAlBKJ1+aeAjoBKwBmitqrtE5AzgbVVtJiK1gBFk/nBJAP6rqs/nsj/r9I0x5jjZfPrGGBNHbD59Y4wxubKib4wxccSKvjHGxBEr+sYYE0es6BtjTByxom+MMXHEir4xxsQRK/rGGBNHrOgbY0wcsaJvjDFxxIq+McbEESv6xhgTR6zoG2NMHLGib4wxccSKvjHGxJGgi76IlBaR6SKyXESm/fGErFzW9YnIQhH5ItjjGWOMCV0onX5vYLqqVgNmZC3npjuQCthTUvIh+0OP4519F3+x7+Iv9l0EL5Si3xwYmfV6JNAip5VE5EzgWmAYkO+nu8Qz+wv9F/su/mLfxV/suwheKEW/nKpuznq9GSiXy3ovAw8DXgjHMsYYUwASj/WhiEwHTs/ho8ezL6iqishRQzcich2wRVUXikhSKEGNMcaELugHo4vIUiBJVTeJSHngK1U9/4h1ngFuBzKAYsCpwFhVvSOH/dl4vzHGBOF4HoweStF/DtiuqgNFpDdQUlVzPZkrIg2Bnqp6fVAHNMYYE7JQxvQHAFeLyHLgqqxlROQMEZmYyzbWzRtjjENBd/rGGGOij/M7ckUkWUSWisgKEXnEdR5XRKSiiHwlIktE5GcRud91Jtfspr5MIlJSRD4RkTQRSRWReq4zuSIij2b9G1ksIqNFpKjrTOEiIu+KyGYRWZztvXzfJPsHp0VfRHzA60AycAFwq4hUd5nJoXSgh6rWAOoB98Xxd/EHu6kv06vAJFWtDtQC0hzncUJEqgB3AXVUtSbgA9q4zBRmw8msldkdz02ygPtO/zJgpaquUdV04APgBseZnFDVTar6Y9brvWT+wz7DbSp37Ka+TCJSAmigqu8CqGqGqu52HMuV38lsjoqLSCJQHNjoNlL4qOosYOcRb+frJtnsXBf9CsD6bMsbst6La1kdTW1grtskTtlNfZnOAraKyHARWSAib4tIcdehXFDVHcCLwDrgV2CXqv7PbSrn8nuT7J9cF/14/7X9KCJyMvAJ0D2r44872W/qI467/CyJQB1giKrWAfaRj1/hY5GInAM8AFQh87fgk0XkNqehIohmXpWTZ011XfQ3AhWzLVcks9uPSyJSBBgLvK+q41zncag+0FxEfgHGAFeJyHuOM7myAdigqt9nLX9C5g+BeHQJMEdVt6tqBvApmX9X4tlmETkdIOsm2S15beC66M8HqopIFRE5AbgFGO84kxMiIsA7QKqqvuI6j0uq+piqVlTVs8g8UfdlTndxxwNV3QSsF5FqWW81BpY4jOTSUqCeiJyY9e+lMZkn+uPZeKBd1ut2QJ7N4jHn3ilsqpohIl2BqWSeiX9HVePyygTgCuDfwE8isjDrvUdVdYrDTJEi3ocBuwGjshqjVUAHx3mcUNVFWb/xzSfzXM8C4C23qcJHRMYADYEyIrIe6EvmTbEfiUgnYA3QOs/92M1ZxhgTP1wP7xhjjAkjK/rGGBNHrOgbY0wcsaJvjDFxxIq+McbEESv6xhgTR6zoG2NMHLGib4wxceT/AZhppPNX9o1QAAAAAElFTkSuQmCC) ### 积分到无穷 In [24]: ```py from numpy import inf interval = [0., inf] def g(x): return np.exp(-x ** 1/2) ``` In [25]: ```py value, max_err = quad(g, *interval) x = np.linspace(0, 10, 50) fig = plt.figure(figsize=(10,3)) p = plt.plot(x, g(x), 'k-') p = plt.fill_between(x, g(x)) plt.annotate(r"$\int_0^{\infty}e^{-x^1/2}dx = $" + "{}".format(value), (4, 0.6), fontsize=16) print "upper bound on error: {:.1e}".format(max_err) ``` ```py upper bound on error: 7.2e-11 ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAlQAAADICAYAAAAuo384AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HX504ISYAkhLAIYScgEBBQFsWFXahVwaWi8KX2URXrUtf+FLWKfrVKq9alrV8rgnX51gXrF4oKohBRwAqigIAQEWQTFNkNCpl7fn/MEBK2JDDJzSTv5+NxHzP3zpl7P3EehnfOOXOuOecQERERkWPnBV2AiIiISLxToBIRERE5TgpUIiIiIsdJgUpERETkOClQiYiIiBwnBSoRERGR41RioDKzCWa22cyWHKXNE2aWZ2aLzKxrbEsUERERqdxK00M1ERh8pBfN7GdAG+dcNnAV8FSMahMRERGJCyUGKufcB8C2ozQ5D/hHtO1/gHQzaxib8kREREQqv1jMoWoCrCuyvx7IisF5RUREROJCQozOYwftH3I/GzPTPW5EREQkbjjnDs43RxSLHqoNQNMi+1nRY4e48cabcM5pi8PtnnvuCbwGbfrsquOmzy9+N3128b2VVSwC1RRgFICZ9QK2O+c2H67h448/xoIFC2JwSREREZHKo8QhPzP7J3AWkGlm64B7gBoAzrmnnXNvmdnPzOxL4AfgV0c+Vx8GDDibLVs2k5AQq9FGERERkWCVmGqcc5eWos11pbmY77/Frl31ueiiX/B///ev0rxFKok+ffoEXYIcI3128U2fX/zSZ1e92LGMEx7ThcxcZK76TGAA//rX6wwbNqxCri0iIiJSFmaGK8Ok9AACFcDlJCa+wnffbSY1NbVCri8iIiJSWmUNVAHdy28CBQV16dt3QDCXFxEREYmhgAKVh++/x8KFn/DYY48FU4KIiIhIjAQ05LffvXjef/PVV6to3rx5hdQhIiIiUpI4mUN1gOd1okmT3axdu7pC6hAREREpSZzMoTrA92exfv1Grr/+t0GXIiIiInJMAu+hinges8uZN28ePXv2rJB6RERERI4k7ob8Drw+kDp1PuH777/VKuoiIiISqLgb8tvPuTfZvbuAoUMvDLoUERERkTKpNIEKEvH9qbz55r95/fXXgy5GREREpNQqzZDfAVdQo8aLrF+/lgYNGpR7XSIiIiIHi9s5VAf4eF42TZo41qz5Es+rRJ1oIiIiUi3E7RyqAzx8/z+sX/8Nw4dfFnQxIiIiIiWqhIEKIBPnpvLaa68yfvz4oIsREREROapKOORX1B143h/5/PMltG/fvlzqEhERETlYFZhDVZznnUrt2iv47rtNJCYmlkNlIiIiIsVVgTlUxfn+LHbvhjPP7Bt0KSIiIiKHVekDFSTh+x/w8cf/4c477wy6GBEREZFDVPohvwPGY3YVM2bMoH///jGrS0RERORgVW4OVXG/IDHx32zYsI7MzMyY1CUiIiJysCoeqHw8rzVNm3p89VWeFv0UERGRclHlJqUXF1n0c+3aDVx22cigixEREREB4i5QATTAucm88srLTJw4MehiREREROJtyK+o2/C8R1i2bCnt2rWL4XlFRESkuov5kJ+ZDTazL8wsz8xuO8zrmWY2zcw+M7PPzezyMtZ8jMYBJ9Oz52ns3bu3Yi4pIiIichhH7aEysxCwAhgAbADmA5c655YXaTMWqOmcG2NmmdH2DZ1zBQedK8Y9VAA/4nmN6dSpJZ999kmMzy0iIiLVVax7qHoAXzrn1jjn9gEvA+cf1OYbIDX6PBX4/uAwVX6S8P1PWLz4cy666BcVc0kRERGRg5QUqJoA64rsr48eK+oZoKOZbQQWATfErrzSaIlz03n99dcZO/beir20iIiICJBQwuulGaO7A/jMOdfHzFoDM8zsJOfcrkObji3yvE90i4U+wFPce+/VtG9/IpdcckmMzisiIiLVQW5uLrm5ucf8/pLmUPUCxjrnBkf3xwC+c25ckTZvAQ845+ZE998DbnPOLTjoXOUwh+pgN+N5T/DRR/Po3r17OV9LREREqqpYz6FaAGSbWQszSwQuAaYc1OYLIpPWMbOGQDvgq9KXHEuP4txAzjjjLDZu3BhMCSIiIlLtlLgOlZkNAR4DQsCzzrkHzWw0gHPu6eg3+yYCzYgEtAedc/97mPNUQA8VRG5P04H09K1s2LCWpKSkCrimiIiIVCVV/F5+pZWP5zWlXbsmfP75Z7rnn4iIiJRJFb+XX2ml4Puf8sUXK7nggouCLkZERESquCoaqACa4dy7TJ48mTFj7gi6GBEREanCquiQX1H/AH7FCy88z8iRIwO4vkj8mzFjBrNmzWLbtm1ce+215OTksHXrVp577jkyMjJo0qQJAwcODLpMEZGY0Ryqw7oNs4eZM+dDTj311IBqEIkfv/3tb3n++ed58MEHOffcc1m0aBHnnHMOq1at4txzz+W///u/mT9/Pvfffz8JCQnk5eWxe/duunbtGnTpIiIxUdZAVdLCnlXEOOALzjqrL1999SVZWVlBFyRSaU2ePJk33niDiRMnMnLkSLKysjj33HMBaN26NZMmTeLkk09mzpw5JCREfoVkZ2czffr0IMsWEQlUFZ5DVZxzbxAOt6Fjx85s37496HJEKq1x48bx61//GoD69evTokUL1q07cAeqRYsWccstt3DZZZexadMmAL777rvCcCUiUh1VkyG//fbieW2oW/cn1q5dTUpKSsD1iFQuq1atIjs7m48//phTTjml8Pibb77JN998w/bt2+ncuTODBg3inXfe4ZFHHqFNmza0bt2am2++OcDKRURiS3OoSpSP57WiUaNEVq/+ksTExKALEqk0nnjiCe666y527NiBWal/j4iIVDlah6pEKfj+F2zatJsTT8yhoKAg6IJEKo3c3Fx69uypMCUiUkbVMFABpOP7y/j660106XIyvu8HXZBIpfDhhx8WG+oTEZHSqaaBCqARvr+EZcvyOPXU04MuRiRwK1euZMuWLXTp0iXoUkq0ZcsWxo8fz1//+tegSxERAap1oAJojnOfMn/+Qvr3HxR0MSKBmjt3LgAnnXRSwJWULDMzk+zsbA3Zi0ilUc0DFUA7nPuIWbNyGTbswqCLEQnM3LlzSUpKom3btkGXElMLFy4MugQRqQYUqADognOzmDx5MqNGXR50MSKBmDdvHu3atcPzqtavhXfeeafw+YQJE3j22WcZNmwYixYtCrAqEalqtBJfod449zYvvDCYOnXq8Ne/Phl0QSIVZteuXSxbtozhw4cHXUqplWbJl82bN9OoUSMA3n77bbp3706nTp3IzMxk1KhRClUiEjNV60/R4zYQeJW//e1vjBlzR9DFiFSYBQsW4Jyjffv2QZdSKlu3bmXq1KnMmTOH1atXH7HdlClTGDZsGAB5eXk8/fTTALRp04Y1a9ZURKkiUk2oh+oQFwITeeihy0lNTWXMmNuDLkik3H388ccAnHjiiQFXUtzrr7+OmbFw4UJycnJ45513mDBhAhkZGTz88MMlvn/r1q2kpaUBcM0117B7924A5syZw5AhQ8q19uMxadIkXnzxRRYuXMiWLVto1qwZF1xwAXfccQe1a9cu8f3r1q3jpptu4t1338U5x4ABA3jsscdo2rRpBVQvUj2ph+qwRgF/4Y477uDBBx8KuhiRcjd//nwAOnbsGHAlByxfvpwzzjiDQYMGMXv2bM4777wyDUmuXr2aVq1aFe4nJCSQnp7O9u3befXVV3nyyZKH9Tdu3MiIESPwPI/Fixcf089xLB555BFq1KjBQw89xLRp0/jNb37DU089xcCBA0sc6szPz6dfv36sXLmS559/nhdeeIG8vDz69u1Lfn5+Bf0EItWPeqiO6BoA7rjjOnbu3MmDD/4h4HpEys+CBQuoUaNGhX7Dz/d9nnjiCcLh8CGvtW7dmqFDhwKRuU/9+vUjJSWFQYNKv7zJm2++yZVXXlnsWDgc5v777+eFF16gfv36JZ6jcePG3HzzzUybNo3OnTuX+trHa+rUqdSrV69w/8wzzyQjI4Nf/vKX5Obm0rdv3yO+95lnnmH16tWsXLmyMFB27tyZ7Oxsnn76aW666aZyr1+kOlKgOqprgFQeeuiX7Nixg7/9TYsIStWzdetW1q5dS4cOHQiFQhV2Xc/zuPHGG4/4+pIlS0hKSuLdd9/l7LPPBuDdd99lwIABpTr/nj17qFmzZrFj//M//8Ott95Ko0aNeOmllxgxYkSJ55k9ezZnnXVWqa4ZK0XD1H77V7DfuHHjUd87ZcoUTj311GK9cy1atKB3795MnjxZgUqknGjIr0QjgX/x1FP/w4gRI4MuRiTmPvvsMwBycnICrqS4t99+m7fffptmzZqxZMkSXnzxRXr27Fmq9y5ZsoROnToVO/baa69x++2306lTJ+rXr8+LL75YqnPNnj2bPn36lLX8mHv//fcBSvziwNKlSw/7WXbo0IFly5aVS20ioh6qUjofmME//3k227fv4M03/x10QSIxs3/hy/IKVCtWrOCFF14gKyuLLVu20KBBA6666qoS3/f//t//O+ZrzpgxgxtuuKHYsYsvvpiLL764xPfOnz+fiRMn0r59e3zfZ86cOdx3332Fr+/Zs4cnn3ySpKQk5s+fz9VXX81HH33EvHnzuO++++jQocMx130kGzZs4O6772bgwIF069btqG23bdtG3bp1DzmekZHBtm3bYl6biEQoUJVaP5yby9tvn84ZZ/Th/fdnVrkFEKV62t9DVR5zhD755BNuvvlm3nrrLWrVqsWqVavK/f57zjnC4fAxDV/Onz+fSy+9lLlz59KgQQMmTJiA7/vFeruefPJJrr/+epKTkxk6dChPP/00EyZM4L777mP06NHFAlVBQQHXXHMN+/btK/Haw4cPLxzaLGr37t2cf/75JCYmMnHixDL/TCJSMRSoyqQ7zi1kzpxT6NatBwsXfqxQJXFv8eLFmBldu3aN+bl/9atf0bt3b1566SV27dpFw4YN+dOf/hTz6xS1cOFCTj/92G54fsUVV3DVVVfRoEEDINLbU3S4zzlH7969SU5OBiK9b48++igJCQns2LHjkPMlJCTw97///ZhqgUhv2LnnnsuaNWt4//33ady4cYnvqVu37mF7orZu3UpGRsYx1yIiR2clfQXXzAYDjwEhYLxzbtxh2vQB/gzUALY45/ocpo2Dklc2jg+r8bxOtGrVlKVLF5GYmBh0QSLHZN++fdSqVYu0tDS+++67mJ575cqVnHjiiWzatKkwoFRmCxYsoEePHixfvpx27doBcM455zBkyBCuu+66Q9pv2LCBli1bsm3bNmrVqhXzevbt28fQoUP58MMPmTFjBj169CjV+/r378/evXv54IMPih3v06cPZsasWbNiXqtIVWRmOOestO2P2kNlZiHgL8AAYAMw38ymOOeWF2mTDvwVONs5t97MMo+t9HjSEt9fyVdfdaRVq7asXLmMlJSUoIsSKbMVK1ZQUFDAySefHPNzb9++HeCQMLV69WpatmwZ8+sdr1WrVpGWllYYpgoKCvjwww8ZN24cH374YWGvl+/7eJ7He++9x8knn1wYpubMmUPv3r2LnXPfvn1ce+21ZR7y832fESNGkJuby9SpU0sdpgDOO+88br311mL/ndesWcPcuXMZN+6Qv4dFJFacc0fcgFOBaUX2bwduP6jNNcB9RztPtJ0DV8W2753nNXCZmY3ctm3bnEi8efXVV52ZuTFjxsT83Hv27HH169d3eXl5hccWLFjgxo0bF/NrxcKSJUtcvXr1Cvcfe+wxV6tWLeecc3/4wx+cc8699tprrmHDhs4554YNG+ZGjRrlnHNu165d7o9//GPMarn66qudmbm77rrLzZs3r9i2fv36wna5ubkuFAq5559/vvDYDz/84Nq0aeM6derkJk+e7CZPnuw6d+7sWrdu7X744YeY1ShS1UUi0tGzTdGtpDlUTYB1RfbXAwd/bzkbqGFms4A6wOPOuReOL+bFiwx8fxVbt3agefNWrFixrPBGrCLxYPnySGdzecyfSkpK4rXXXuP++++nV69e7N27lyZNmhzXt/fKU05ODjfffDP33XcfqampdO3alb59+/LHP/6R0047DYCsrCzOPPNMHnnkEW655RaefPJJnnrqKfLz87n++utjVsu0adMwMx544AEeeOCBYq+NHTuWu+++Gyj+B/F+KSkpzJw5k5tuuon/+q//KnbrGfWki5Sfo86hMrMLgcHOuSuj+yOBns6564u0+QvQDegPpADzgHOcc3kHncvBPUWO9IluVcFePK8ziYnrmDPngxK/1ixSWVx66aW88sor5OXl0bp166DLEREJTG5uLrm5uYX79957b5nmUJUUqHoBY51zg6P7YwDfFZmYbma3AcnOubHR/fFEhgknHXSuKjQp/XB8PG8wMJNXXnmZiy66KOiCRErUpUsX1q1bx/fffx90KSIilUpZJ6WX9J3/BUC2mbUws0TgEmDKQW0mA6ebWcjMUogMCVbD5Xg9fP8dfP9aLr74F8UWAhSpjHzfZ8WKFaVefVxERI7sqHOonHMFZnYdMJ3IsgnPOueWm9no6OtPO+e+MLNpwGLAB55xzlXDQLXf40AH7rnnGpYsWcZrr70cdEEih7V69Wp++uknevXqFXQpIiJxr8R1qGJ2oSo/5HewXMzOJienIwsWfKS1qqTSmTx5MsOGDWP69OkMHDgw6HJERCqVWA/5yTHrg3NfsHTpGpo0ac6mTZuCLkikmKVLl2JmGvITEYkBBapy1RLfX8u2bXVo0aJV4U1oRSqDRYsWkZOTQ2pqatCliIjEPQWqclebcPgL9u49g+7de/Dqq68GXZAIAJ9++ilnnXVW0GWIiFQJClQVwsO56fj+9VxyyXDGjr036IKkmtu+fTurVq3izDPPDLoUEZEqQYGqQv0Z+Dv33nsfw4ZdiO/7QRck1dTs2bMxM/r37x90KSIiVYK+5ReI2ZidTVZWIxYs+M8hN48VKW+jR49m5cqVzJo1q/DYQw89RNu2bVm4cCGjRo2ibdu2AVYoIhIsfcsvLpyJc+vYsMEjK6sZ06dPD7ogqeJ83+fss8/mgw8+YPfu3UyaNInRo0cXvj5nzhxWrlzJBRdcwG9+8xt+97vfBVitiEj8UaAKTCa+n8e+fRcxePAQbr1V/4BJ+dm9eze5ubls3LiR22+/ndatWzN8+PDC12fNmkWPHj0AaNKkCfPnzw+qVBGRuKRAFSgPeBF4jkcf/TNdu55Cfn5+0EVJFZSamsoDDzzAVVddxcKFC5k0qditNtm8eTMpKSmF+6FQiO3bt1d0mSIicUuBqlIYhXPLWbz4axo2PIHPPvss6IKkCrr11lvZsWMHc+fOpVmzZsVe832fUChUuF9QUFBsX0REjk6BqtLIxve/IT//ZLp1O5nHH3886IKkGmnSpAk//PBD4X44HKZOnToBViQiEl8UqCqVBHx/Js7dx4033syQIedoaQWpEAMGDCjsGc3Ly6N79+4BVyQiEl+0bEKlNQfPO5vMzFTmz//okCEakVi744476NSpE59++ilXXnkl2dnZQZckIhKYsi6boEBVqe3E807DbCXPPTeBkSNHBl2QiIhItaBAVSX9FvgLffv2Z9q0N0lMTAy6IBERkSpNgarK+g+eN4SkpALefHMKffr0CbogERGRKksrpVdZPfH9b9mzpx99+/bj8st/pQnrIiIilYR6qOLSG5hdRv36dcnNfY/27dsHXZCIiEiVoh6qamEYzm1my5amdOyYwz33jA26IBERkWpNPVRx73HMbqFt23bMnj2LBg0aBF2QiIhI3FMPVbVzA859RV7ejzRunMWzzz4bdEEiIiLVjnqoqpRbgD/Tq9dpvP32VNLT04MuSEREJC6ph6paewRYwPz5a8jMbMCjjz4adEEiIiLVgnqoqqw7MRtHixatmD79Td1GREREpAzUQyVRD+DcWr7+OoV27U5k9OjfaN0qERGRclJioDKzwWb2hZnlmdltR2nX3cwKzOyC2JYox64xvv8Zzk1g/PjnyMioz8yZM4MuSkREpMo5aqAysxDwF2Aw0AG41MwOWUUy2m4cMA0odfeYVJRf4vvb2LnzVPr3H8CgQUPIz88PuigREZEqo6Qeqh7Al865Nc65fcDLwPmHaXc9MAn4Lsb1Scwk4dxUIJf33ltA3br1tMSCiIhIjJQUqJoA64rsr48eK2RmTYiErKeihzTzvFI7E9/fzN69V3DFFVfRoUNnVq1aFXRRIiIicS2hhNdLE44eA253zjkzM4465De2yPM+0U0qngc8CdzIihXnk52dzdChF/Dii8+TkpISdHEiIiIVLjc3l9zc3GN+/1GXTTCzXsBY59zg6P4YwHfOjSvS5isOhKhMIB+40jk35aBzadmESusNPO8KPG83v//9ndx9991BFyQiIhKosi6bUFKgSgBWAP2BjcDHwKXOueVHaD8R+Ldz7l+HeU2BqlLzgfsw+wNpaak899yznH/+4abLiYiIVH0xXYfKOVcAXAdMB5YBrzjnlpvZaDMbfXylSuXiAWNxbivbt5/J0KHDaN++EytWrAi6MBERkUpPK6XLEeTheRfi+59z7rnn8b//+yK1a9cOuigREZEKoZXSJUay8f3FwP/x1ltzSU/P4M4779Rq6yIiIoehQCUlOI9w+FvC4bt48MGHSU2ty+OPPx50USIiIpWKhvykDH4ErsfsOVJTU3n44XFcccUVQRclIiIScxryk3KUBDyDczvYseMcrrrqaurVa8hLL70UdGEiIiKBUqCSY5ACPI9zW9i69QxGjhxFw4ZNeOONN4IuTEREJBAKVHIc0oncwnEz333XlQsuuJCsrBZMnz496MJEREQqlAKVxEBm9MbL69m4MZvBg4fQsmU2s2fPDrowERGRCqFAJTHUGOdmAKv5+usTOOusPjRt2oJXX3016MJERETKlQKVlIPmODcbWMOGDR245JJLqVs3k0cffVTrWImISJWkZROkAuwkstzCP6lZM5HrrvsNDzzwAImJiUEXJiIiclgxvTlyLClQCewFfo/n/RWzvQwfPpy//OUJ0tPTgy5MRESkGK1DJZVYIjAO399JOPwnXn55GhkZ9Rg0aDBr164NujgREZFjpkAlAfCAGwiHv8W5V5g5cwXNm7cgJ+ck3nrrraCLExERKTMFKgnYRYTDq4E5LFtWm3PO+TlpafUYM2YMP/74Y9DFiYiIlIrmUEklsxO4Hc97AdhD3779ePzxP9OxY8egCxMRkWpEc6gkzqUCf8P3d+H7/yA3dw05OZ1o1qwV48eP17ILIiJSKamHSuLAKuAGzKaTmFiDyy4bzsMPP0xGRkbQhYmISBWlHiqpgloDU3FuDz/9dBvPPz+VevUyadeug3qtRESkUlAPlcSp/2D2eyCXUMgYMKA/Dz74B7p06RJ0YSIiUgWoh0qqiZ449w7O/UhBwZ+ZMeNLunbtRr16Dfnd737H7t27gy5QRESqEfVQSRWyCfg9odBrhMM7ycnpzF13jeGSSy4JujAREYkz6qGSaqwR8Azh8HZgJkuXpnHppSOoWTOZwYOHMHv27KALFBGRKko9VFLFFQBP4Hnj8f0vSEpKYdCgAfz+93dxyimnBF2ciIhUUro5ssgR5QN/JhT6B+Hwl9SqVYdzzhnC3Xf/XguHiohIMeUy5Gdmg83sCzPLM7PbDvP6CDNbZGaLzWyOmXUuS9EiFSMFuJNweCWwnR9++C2vv/4ROTk5pKXVY9SoX7Jq1aqgixQRkThUYg+VmYWAFcAAYAMwH7jUObe8SJtTgWXOuR1mNhgY65zrddB51EMlldQW4EFCoVcIhzeQmprB2Wf355ZbbqFnz55BFyciIgGI+ZBfNCzd45wbHN2/HcA599AR2tcFljjnsg46rkAlcWAT8Aih0L8Ih1dTs2YSvXr14pprruaiiy7C8/Q9DhGR6qA8hvyaAOuK7K+PHjuSXwNvlbYAkcqlEfAnwuFVQD4//fQHPvhgO8OHj6BGjZp06tSFRx55hPz8/KALFRGRSqQ0PVQXAoOdc1dG90cCPZ1z1x+mbV/gr0Bv59y2g15zcE+RI32im0g88IE3MHsKs3n4/h6ysppz0UVDue6662jdunXQBYqIyHHIzc0lNze3cP/ee++N+ZBfLyJzovYP+Y0BfOfcuIPadQb+RSR8fXmY82jIT6qQBUSGBmcSDn9LUlItunc/hZEjL2PUqFEkJSUFXaCIiByH8phDlUBkUnp/YCPwMYdOSm8GzARGOuc+OsJ5FKikitoJPIPZK5gtwfd/olGjJgwa1J/rrruW7t27B12giIiUUbmsQ2VmQ4DHgBDwrHPuQTMbDeCce9rMxgPDgLXRt+xzzvU46BwKVFJNLAT+Rig0g3B4HQkJieTk5PCLX1zIlVdeSWZmZtAFiohICbSwp0ilshd4CbMX8LwFhMO7SEmpQ5cuJzF06Hn8+te/JiMjI+giRUTkIApUIpXat8AEzP6N5y0mHN5NrVqpdOlyEhdcMJTLL79cAUtEpBJQoBKJK5uAZzGbiuctIRz+gVq10ujWrQvnn38uI0aMoFGjRkEXKSJS7ShQicS1jRzowVpOOLyLxMRk2rRpQ79+ZzF8+HBOPfVULTAqIlLOFKhEqpSdwMvAv0lIWEhBwSbMoEGDE+jVqzvnn38eF198MbVr1w66UBGRKkWBSqRK84HZwCt43gfAKnz/R1JSUsnObkPv3r0YOnQoffv2JSEhIeBaRUTilwKVSLWzHvgn8B4JCZ8TDm/COZ86ddI58cS2nHXWGQwdOlRDhSIiZaBAJSLAUmAS8D6h0FLC4S2YOdLSMsnJac9pp/XinHPO4fTTT1fIEhE5DAUqETkMH/iESMj6kISEPMLh73HOp1atVJo3b0737t3o168fP//5z7V0g4hUewpUIlIGK4ApwGxCoWU4txHf/5GEhJo0bNiIzp07cNpppzFw4EC6d++u3iwRqTYUqETkOO0E3gbexfM+wWwN4fAOwCcpqRYnnNCYnJz29OzZg0GDBnHyyScraIlIlaPuEvbiAAAKXElEQVRAJSLlZBWRoDUPz/scs3WFQSs5uXZh0OrWrSunn346vXv3JikpKeCaRUSOjQKViFSwPGA6MBfP+xzP20g4vB3nwiQkJJKeXo8WLZqSk9OBHj160K9fP7Kzs9WrJSKVmgKViFQSW4GZwDxgMaHQV8BmwuF8AJKTa5GZ2YCWLZvSvv2JdO3ald69e9OhQweFLREJnAKViFRyPrAcmAUsBPJISFiH72/B9/MBR40aSaSn16Vp0ya0a5dNp06d6NmzJz169NCq8CJSIRSoRCTOrQU+ILLMwzJCoTXAZnx/F86FMfNITq5F3boZNG3ahNatW9KxY0e6dOlCz549teSDiMSEApWIVGE/EglaC4ElwCpCoXXAlmjgKsDMo2bNZNLS0mnYsAEtWjSlVatWtG/fnpNOOomTTjpJk+VFpEQKVCJSjRUAn7G/dwu+wvPW43nf4dx2fH8PzvnR0JVCWloajRo1pGnTxrRo0YI2bdpw4oknctJJJ9GoUaNgfxQRCZQClYjIUeUDnwKLiczlWoXnfYPnbcG5Hfh+Ps4VABAK1SA5uRZpaWk0bFifrKzGNG/evDB8tW/fntatW2sSvUgVpEAlInLcfGAdkWHFSOiCr/G8TdHgtQvn9uD7+wCHmUeNGjVJSalNenoaDRpk0rhxI7KysmjWrBktW7akTZs2tG3blpSUlCB/MBEpJQUqEZEKtZ1I6FpJJHitBTbieZvxvG0HhS8fMEKhBBITk0hJqUV6ehqZmRk0atSARo0accIJJ5CVlUWLFi1o1aoVTZs2JSEhIcCfT6R6UqASEam0CoCviYSv1dHnG4DNeN4WPG8b8EM0gP1UOPR4IITVJDk5hTp16lC3bjoZGenUq5dB/fr1OeGEEwrDWLNmzWjevLl6w0SOgwKViEiVkk8kfK0m0vu1HvgG+BbYRii0A7NdQD7O/Yjv740GscjvW88LEQrVIDGxJklJydSqlUKdOrWpWzeN9PQ06tWrR7169WjYsCENGzakcePGNG7cmCZNmpCamqr5YVJtKVCJiAiR4cWNwBoi88E2R7ctwPfAdjwvEsbMIr1izv2Ec/uKBTKwaChLoEaNRGrWPBDMUlNrk5pah7S0VNLS0khPT6du3brUq1ePzMxMMjMzadiwIY0aNSIzM1PhTOKKApWIiMTITiK9Yd8QCWPfAt8RCWTbiMwf24nn7cYsH7M9wE9Fgtk+nPOJhLsIM68woCUkJJCQsD+kJZGcnETt2inUrl2L2rVrkZqaSu3atalTpw5paWmFW0ZGBnXr1iUjI4OMjAwyMzNJTEys6P84UsUpUImISCW0l0gg20QkkG0lEsq2AjuIhLNdRELcD5j9gOflY/YT8BOwF+f2AvtwLoxzBdGwVvTfFcPM8LxIaPO8/aGtBomJNUhMjIS35OSkwi0lJZnk5GRSUlIKt9q1axcGuTp16hQGu/2P+3vjUlJS1OtWhcU8UJnZYOAxIASMd86NO0ybJ4AhRAb7L3fOfXqYNgpUcS0X6BNwDXJsctFnF89y0edXEp9IIPs++ridSEjbv+0Cdke3H6KP+cCPmO3ffsJsL5HgtxcoiIa2AiCMc2HAJ/Jv5qH/lpl5mFlhD5zneTjnSEysWaQ3LoEaNYqGu0QSE2uQlFQz+vxAb93+x0jPXXLhY9HnSUlJpKSkkJycTK1atQoDYdHnCnzHrqyB6qjfxTWzEPAXYACRr6LMN7MpzrnlRdr8DGjjnMs2s57AU0CvY6peKrFc9Es9XuWizy6e5aLPryQekBHdysa5yHZsCoj0qO2ILo+xE9hJOLw/uE1i794+RMJbJMDBniKPe6PP9wL5mG3D8/YB+zArYH+wgzD7A14kPIYLh1IP9NIdOewdEOnB2/94IAB6hY8HthChkEcoFBmeDYU8EhISCIVC0XAYokaNGtHH/WGxeHAsuoVCocLnCQkJJCYmRsNlYuF+0ef72xZ9XrNmzWLPi55j/2PRLSEhoUIDZUmLm/QAvnTOrQEws5eB84ksurLfecA/AJxz/zGzdDNr6JzbXA71ioiIVBIJHD3IfQ38rtRncw7C4RiUVcgnEtj2B7o9OLenyGPRgPcj+4dWiz8vuu076LEg+rygyPYjZgWYhaOhMBIGzfzo8zD7Q2HksXhAjGwHB8TIdvj9/cqSiiOdTmZFn+8PmnAgeJZNSYGqCZGvh+y3HuhZijZZRGYwFpOaem6ZC5TK4ccfV5CU9EnQZcgx0GcX3/T5xa/q8dmFolvNoAs5RCSkFUSHaw/08EWGciPhLjIfL1xsH3zC4e+JhMbSKylQlTbyHRzlDvu+nTunlvJ0Uhnt3ZsXdAlyjPTZxTd9fvFLn131UVKg2gA0LbLflEgP1NHaZEWPFVOWiV0iIiIi8aSk2VoLgGwza2FmicAlwJSD2kwBRgGYWS9gu+ZPiYiISHVy1B4q51yBmV0HTCcySPqsc265mY2Ovv60c+4tM/uZmX1J5GsNvyr3qkVEREQqkQpb2FNERESkqir3BRrMbLCZfWFmeWZ2W3lfT2LHzJqa2SwzW2pmn5vZb4OuScrOzEJm9qmZ/TvoWqT0okvQTDKz5Wa2LDqlQuKEmY2J/u5cYmb/a2aV72twUsjMJpjZZjNbUuRYhpnNMLOVZvaOmaUf7RzlGqiKLAw6GOgAXGpm7cvzmhJT+4CbnHMdiSzWeq0+v7h0A7AM3aog3jwOvOWcaw90pvj6f1KJmVkL4Eqgm3OuE5EpM8ODrElKNJFIVinqdmCGc64t8F50/4jKu4eqcGFQ59w+YP/CoBIHnHObnHOfRZ/vJvILvXGwVUlZmFkW8DNgPIcubyKVlJmlAWc45yZAZD6rc25HwGVJ6e0k8gdpipklACkc5tvvUnk45z4gcnPJogoXLo8+Dj3aOco7UB1u0c8m5XxNKQfRv7i6Av8JthIpoz8TWarZD7oQKZOWwHdmNtHMFprZM2aWEnRRUjrOua3AI8BaYCORb7+/G2xVcgyK3vVlM9DwaI3LO1BpiKEKMLPawCTghmhPlcQBM/s58G30ZuXqnYovCUA34G/OuW5EvkF91OEGqTzMrDVwI9CCSK9+bTMbEWhRclxcyTdKLPdAVZqFQaUSM7MawOvAi865/wu6HimT04DzzGw18E+gn5k9H3BNUjrrgfXOufnR/UlEApbEh1OAuc65713kfif/IvL/o8SXzWbWCMDMTgC+PVrj8g5UpVkYVCopi9wd8llgmXPusaDrkbJxzt3hnGvqnGtJZELsTOfcqKDrkpI55zYB68ysbfTQAGBpgCVJ2XwB9DKz5Ojv0QFEvhgi8WUK8Mvo818CR+1UKOnWM8flSAuDluc1JaZ6AyOBxWb2afTYGOfctABrkmOnIfj4cj3wUvSP0VVo0eS44ZxbFO0NXkBk/uJC4O/BViVHY2b/BM4CMs1sHXA38BDwqpn9GlgD/OKo59DCniIiIiLHp9wX9hQRERGp6hSoRERERI6TApWIiIjIcVKgEhERETlOClQiIiIix0mBSkREROQ4KVCJiIiIHKf/D3hzxzX6sL0TAAAAAElFTkSuQmCC) ### 双重积分 假设我们要进行如下的积分: $$ I_n = \int \limits_0^{\infty} \int \limits_1^{\infty} \frac{e^{-xt}}{t^n}dt dx = \frac{1}{n}$$In [26]: ```py def h(x, t, n): """core function, takes x, t, n""" return np.exp(-x * t) / (t ** n) ``` 一种方式是调用两次 `quad` 函数,不过这里 `quad` 的返回值不能向量化,所以使用了修饰符 `vectorize` 将其向量化: In [27]: ```py from numpy import vectorize @vectorize def int_h_dx(t, n): """Time integrand of h(x).""" return quad(h, 0, np.inf, args=(t, n))[0] ``` In [28]: ```py @vectorize def I_n(n): return quad(int_h_dx, 1, np.inf, args=(n)) ``` In [29]: ```py I_n([0.5, 1.0, 2.0, 5]) ``` Out[29]: ```py (array([ 1.97, 1\. , 0.5 , 0.2 ]), array([ 9.804e-13, 1.110e-14, 5.551e-15, 2.220e-15])) ``` 或者直接调用 `dblquad` 函数,并将积分参数传入,传入方式有多种,后传入的先进行积分: In [30]: ```py from scipy.integrate import dblquad @vectorize def I(n): """Same as I_n, but using the built-in dblquad""" x_lower = 0 x_upper = np.inf return dblquad(h, lambda t_lower: 1, lambda t_upper: np.inf, x_lower, x_upper, args=(n,)) ``` In [31]: ```py I_n([0.5, 1.0, 2.0, 5]) ``` Out[31]: ```py (array([ 1.97, 1\. , 0.5 , 0.2 ]), array([ 9.804e-13, 1.110e-14, 5.551e-15, 2.220e-15])) ``` ## 采样点积分 ### trapz 方法 和 simps 方法 In [32]: ```py from scipy.integrate import trapz, simps ``` `sin` 函数, `100` 个采样点和 `5` 个采样点: In [33]: ```py x_s = np.linspace(0, np.pi, 5) y_s = np.sin(x_s) x = np.linspace(0, np.pi, 100) y = np.sin(x) ``` In [34]: ```py p = plt.plot(x, y, 'k:') p = plt.plot(x_s, y_s, 'k+-') p = plt.fill_between(x_s, y_s, color="gray") ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVfXa//H3VwbJEXFONCewHHFi2iioleijZc5alkpP46mec9kvw9JMq2OWdTyp2DH1yk6KRT5qZWIqVGqS5aygR0RFHFLxoCDz/v7+EHmMQKa9WXu4X9fldfaGtdf+tDzeLO691vdWWmuEEEI4llpGBxBCCGF5UtyFEMIBSXEXQggHJMVdCCEckBR3IYRwQFLchRDCAZVb3JVSK5RSF5VSh+6wzT+UUv9WSh1QSvW0bEQhhBCVVZEz95VAeFnfVEoNBTpqrX2Ap4AoC2UTQghRReUWd631T8DVO2zyEPBp0bYJgKdSqrll4gkhhKgKS/TcWwGptz0/C3hbYL9CCCGqyFIfqKoSz2VNAyGEMJCrBfaRBrS+7bl30df+QCklBV8IIapAa13yBLpcljhz3wg8DqCUCgT+o7W+WNqGWmu7/fPGG28YnsFR81+6dIlZs2b94XlcXFyp2dPT0+nRowcdO3bk1VdfJTQ0lOHDh1OvXj1iYmJK3f/Jkyc5fPhw8fPly5cTHx8vx17y28Wfqir3zF0ptQYIBZoopVKBNwC3omL9sdZ6k1JqqFLqBJAFTKlyGuE0zpw5Q+vWrVFK0bBhQ1q2bInWGqUUTZo0ISws7E+vOXHiBAMHDqRp06ZMmDABFxcXAHr37k3Dhg154oknSElJ4eWXX/7D69q1a/eH576+vjRr1qz4eWJiIp06daJWLbntQziOilwtM0FrfbfW2l1r3VprvaKoqH982zZ/0Vp31Fr30FrvtW5k4QgmT57MiRMnAHBzc+OZZ55BqbJ/8/zpp5/w9/fHx8eH4cOHFxf2tm3bAtCxY0cmTZrE3LlzefbZZzGbzWXuKyQkBF9fX+Dmb5MvvvgiaWl/6iQKYdfkVKWCSjuTtCdG59+3bx8//vhj8fNt27bh4+NTodcWFBQwZMgQwsLC6N+//x9+CNx+Vt6iRQsiIiJYv349Q4cOJS8vr9x9K6X4/vvvad365sdGFy5cYOnSpRX9z6oQo499dUl++6Sq09Op1BsppWvqvYTt2b59OxkZGTzyyCOVet1bb73FvHnzGD169J/aK2XJyckhJiaGOnXqEBcXh5eXV4XfLyUlhbi4OKZOnVqpnEJYi1IKXYUPVKW4C6vIy8vjww8/ZNq0abi6Vv6iLLPZTEREBBs2bGDChAl/6JFXRGFhIZs2beLixYts27atuA1TWQsXLmTIkCFVfr0Q1VXV4i5tGWEVbm5uuLi4kJ2dXenX5uTkcP/997NlyxamTp1a6cIO4OLiwrBhw+jUqRMBAQH88MMPld4HQNOmTSt15i+ErZAzd2ExFy5c4NixY4SGhlZ5HxcvXmTAgAForRk5ciTu7u7VznXw4EFiY2OJiopi0qRJVd7PyZMnSUpKYujQodXOJERFyZm7MFxqaiq7d++u8uuPHDlCz549adCgAWPHjrVIYQfo3r07Y8aM4dlnn2XOnDlV3s+VK1e4eLHUWziEsDly5i6q5caNG7i4uFC7du1q7WfLli2MGTOGoKAgAgMD73hZZFVdunSJ1atXM3z4cFauXFnt69rPnj2Lt7csoySsS87chSFmzpzJF198Ua19LFu2jJEjRxIeHk5QUJBVCjvc7J9HRESwbds2Bg4cWKXPA265cOECo0aNorCw0IIJhbAcOXMX1ZKbm4u7u3uVC3JkZCSLFi1i3LhxxdeaW1teXh7r1q0DIC4ujhYtWlRpP2azWe5qFVYnZ+6ixsyfP5/k5GQAateuXaXCXlhYyLhx41i2bBlTpkypscIO4O7uztixY/H09KRnz54cPny4Svu5Vdjz8vJ48sknuXr1TmMPhKhZUtxFpfn4+FCvXr0qvz4zM5N+/fqRkJDA1KlTady4sQXTVUytWrUYPHgwfn5+mEwmYmNjq7wvNzc3hg8fjqenpwUTClE90pYRFZKdnc1dd91V7f2kpqYycOBAPDw8ePjhh3Fzc7NAuuo5evQo3377LR988AFPPfVUtfeXnp4u18YLi5G2jLCqhx9+mH379lVrH3v37qV37940b96ckSNH2kRhB+jcuTPjx49n2rRpTJ8+vVr7KigoYNCgQXLJpDCcnLmLCsnIyKBhw4ZVfv3XX3/No48+Sv/+/enbt68Fk1lOeno6q1evJiwsjOjo6OKVJysrLy/PYtfoCyFn7sLi9u/fX3y5YHUK+0cffcSECRMYPny4zRZ2AC8vL6ZOncqePXsICQkhMzOzSvu5Vdi11nz99dd3XH5YCGuR4i7KtGLFCg4cOFCtfbz00ku8/vrrPProo3ax+FadOnV47LHHyMrKws/Pj9TU1PJfVIbc3FzWr19f5R8SQlSHtGWEVeTn5zNq1CgSEhKYOHGi3V1JYjab2b59O0lJSWzevJnevXsbHUk4KWnLCIv44Ycfqn22npGRQWBgIEeOHGHKlCl2V9jh5qWS999/P4GBgQwYMIANGzZUa39Xr17l7bfflhaNqDFS3MUfXLp0qVo345w8eZIePXpgNpuZOHGiRS6fNFKfPn0YPnw4EydO5O9//3uV9+Pm5kbjxo2ttrSCECVJW0ZYzK5duxg+fDjdu3cnNDTUoQrZuXPniI6O5vHHH2fhwoWy7ICoMdKWEVWWlJTEkiVLqrWPtWvXMnjwYPr3709YWJhDFXaAu+++m6lTp/Lll1/y0EMPkZ+fX+V97dixgwULFlgwnRB/JsVdULduXZo3b17l18+bN4+IiAhGjhyJn5+fBZPZFk9PT6ZMmUJSUhIBAQFVbl917NiRwMBAC6cT4o+kLSOqzGw289RTTxETE8PEiROr9QPCnhQUFPDtt99y5coVtm/fTocOHYyOJByYtGVEpeTm5jJ9+nSysrKq9PqcnBwGDx7Mpk2biIiIcJrCDuDq6spDDz1Ehw4d6Nu3Lzt37qzyvmbPnk1CQoIF0wlxkxR3J+Xq6kqnTp2qdDXL5cuX6du3L2fOnGHy5Mk0aNDACgltm1KK0NBQQkNDGTx4MNHR0VXaz4MPPmgXN3cJ+yNtGVEpiYmJPPDAA9x9992Eh4dXef0VR3Ly5Em++uorIiMjee2114yOIxyMtGVEhaxdu7bKQ6zj4uIICgqic+fODB06VAp7kfbt2/PEE0/w/vvvM3Xq1CrdqHTlyhWeffbZal2FI8TtpLg7mQYNGlC/fv1Kv27lypUMHz6cBx54AJPJ5HCXOlZXs2bNiIiIIDY2lgceeICcnJxKvb5Ro0Y88MADuLq6WimhcDbSlhHlmjVrFh9++CFjxozhnnvuMTqOTcvNzWXdunW4uLgQHx9P06ZNjY4k7Jy0ZUSZsrOz+fjjj6nsD9dbSwgsXryYyZMnS2GvgNq1azN27Fjq1atHjx49SExMrPQ+Nm7cyLZt26yQTjgTKe5OIDMzk8uXL1fqNVlZWYSGhrJjxw4iIiJo0qSJldI5HhcXF4YMGUK3bt0ICgpi69atlXp9o0aN7HKxNWFbpC0j/uTcuXMMGDAANzc3RowYYTPj8OzRkSNH2LRpE//4xz+YOnWq0XGEHZK2jPiTvXv38vvvv1fqNfv376dnz540adKEUaNGSWGvpi5dujBu3LjioSWVkZuby6efflrpdpoQIMXdoW3fvr1SQ62//fZb+vfvT58+fXjwwQdl5UMLadOmDVOmTGHp0qWMHz+ewsLCCr2uoKCAQ4cOkZeXZ+WEwhFJW0YAsGTJEl555RWGDx/Ovffea3Qch5SVlcXatWvx9vZmy5Yt1K1b1+hIwg5IW0YUO3v2bKW2nzZtGq+++ioTJ06Uwm5FdevWZdKkSWRkZNCzZ0/S0tIq/Nrk5GSuXbtmxXTC0ZRb3JVS4UqpJKXUv5VS00v5fhOl1Gal1H6l1GGl1GSrJBUVcu3aNR566CGys7PL3bawsJARI0bw2WefMXXqVFq1alUDCZ2bm5sbo0aNolmzZvTq1avCbbN//vOf/Pzzz1ZOJxzJHdsySikX4BhwP5AG7AEmaK0Tb9tmNlBbax2plGpStH1zrXVBiX1JW6aGmM3mcvvl165dY9CgQVy5coWxY8fa/Tg8e7Rnzx5+/PFHVq9ezbBhw4yOI2yUtdoy/sAJrfUprXU+EA08XGKb88CtZQEbAFdKFnZhfVrr4jVNyivsp06dws/Pj/z8fB599FEp7Abp27cvw4YNY/z48Xz00UcVft2NGzesmEo4ivKKeysg9bbnZ4u+drtlQBel1DngAPCS5eKJivriiy94+eWXy90uISGBPn360KZNG0aMGCFrmRisU6dOTJw4kddee43/+Z//KXf7H3/8kUcffbQGkgl7V96/7Ir0UWYA+7XWYUqpDsD3SqkeWuvrJTecPXt28eOwsDDCwsIqEVXcyahRoxgwYMAdt/nqq6+YPHkyAwcOpFevXjWUTJSnVatWTJ06ldWrV5OSkkJMTEyZ9xf069eP3r1713BCUZPi4+OJj4+v9n7K67kHArO11uFFzyMBs9b63du22QS8rbXeWfR8GzBda/1riX1Jz91A77//PrNnz+aRRx6hY8eORscRpcjOzubLL7/Ey8uLrVu30rBhQ6MjCRtgrZ77r4CPUqqtUsodGAdsLLFNEjc/cEUp1RzoBJysbBBRNfPnz2fv3r1lft9sNvPMM8/w1ltvMWnSJCnsNuyuu+5i4sSJFBQU4OfnR0pKSpnb5uTkMGLECLk8UpTpjsW96IPRvwCxwFFgrdY6USn1tFLq6aLN3gH6KKUOAFuBV7TW6dYMLf5Pz549adOmTanfy8vLY+jQoWzYsIGIiAhatGhRw+lEZbm6uvLwww9zzz330KdPnzIHq3h4ePD//t//q9La/MI5yB2qDurKlSsMHDiQ7OxsRo0ahYeHh9GRRCXt3buXbdu2sXLlSsaMGWN0HGEQuUPVyWzatKnMNUqOHz+On58fbm5ujB8/Xgq7nerVqxcjR45kypQpzJ8/v8ztYmJiOH/+fA0mE/ZAirsdys3NJTo6utTrnX/44QcCAgLo1KkTw4YNkzmndq5Dhw48/vjjvPPOOzz11FOlzme9dOkSV65cMSCdsGXSlnEgn332Gc8++yyDBw+me/fuRscRFnTt2jXWrFlD165d+eabb3B3dzc6kqghVW3LSHG3M1lZWaWuJjhnzhzmz5/PmDFjaNu2bc0HE1aXm5tLTEwMtWvXJi4u7k/TsXJycrhx4wZeXl4GJRTWID13J3Ds2DEefPDBPwxvMJvNPP744yxcuJDJkydLYXdgtWvXLv4Mxc/Pj6SkpD98f+nSpSxfvtygdMLWyJm7nbn9zD07O5vw8HBOnjzJ+PHjqVevnsHpRE3QWrNz50727NnD+vXri+/0LiwspFatWihV6ZM8YcOkLeNkLly4wIABA1BK8cgjj0gP1gkdOnSIzZs3s2jRIiZPnmx0HGEl0pZxYDt37mTBggXFzw8ePEjPnj3x9PRkzJgxUtidVLdu3Rg7diwvvPDCH9Zt2rVrFy+9JOv3OTtZEtAOtG/fvvhxbGwsY8eOJTg4mMDAQANTCVtwzz33MHnyZD766CNOnDjBqlWr6NGjBw0aNCj/xcKhSVvGjixdupSXX36ZYcOGcd999xkdR9iQzMxM1q5dS9u2bYmNjZU1+h2ItGUcUEpKSvGczenTp/PKK68wYcIEKeziT+rVq8ekSZO4cuUKfn5+nDt3juvXr7NxY8l1/oSzkOJuw3766SfWr1/P6NGjWbFiBVOmTMHb29voWMJGubu7M3r0aBo3bkyvXr349ddf2b59O/Ibs3OStowNy8zM5P777+fixYuMHTuWOnXqGB1J2ImEhAR27txJdHQ0Q4YMMTqOqAZpyziQuLg4zpw5Q48ePbhx4waPPvqoFHZRKQEBAQwZMoQxY8awdOlScnJyLDLdR9gPOXO3Qd26deP06dP4+fkxYMCAcgdeC1GWs2fPEh0djbu7O8899xxz5841OpKoJDlzdxDbtm3j2LFjmEwmBg0aJIVdVIu3tzcREREopVi/fr3RcUQNkjN3G3FrKO7ixYu5fPkyoaGhALRt25Z27doZnE7Yq5SUFE6dOkV+fj67du0iIiICb29vGVBvR6p65i43MdmIsLAwtm/fTm5uLv369WPAgAFGRxIOoF27dsUnB2fPnuX7778nOTkZV1f5p+/o5Hd+G/Lpp5/So0cPacUIq/D29ub8+fPs37/f6CiiBkgVsRE7duwgPT2d/v37y7K9wip8fHzw9/fnzTffNDqKqAHSc7cBN27cYMiQIbi5udGvXz+j4wgHdv36dZYsWUJcXBwBAQFGxxEVIFfL2LERI0aQkJBA3759jY4iHFz9+vXp0qUL4eHhXL161eg4woqkuNsApRR9+/bFw8PD6CjCCQQHB5Obm0teXp7RUYQVSXE3WHJyMj/99JMs3ytqjJeXF76+vsycOdPoKMKKpLgb6OOPP+b555+nR48eMiJP1CiTycSaNWt47LHHuH79utFxhBVIcTdQfn4+P/zwg5y1ixrXrFkzWrduTWZmpkzyclBS3A20d+9eOnfujKenp9FRhBMymUzEx8djNpuNjiKsQIq7Aa5cuUJ6ejpffPEFwcHBRscRTsrb25smTZowf/58jh07ZnQcYWFynbsBHnvsMTIzM0lOTmb06NFGxxFO7OTJk3z33Xd07NiRTZs2yW+RNkjWlrEjixcvpk2bNowfP97oKMLJtWvXDg8PD0aPHi2F3cFIW8YA77//Ps2bN+fuu+82OopwckopTCYTH374ofTeHYwU9xq0fPlydu3aRVRUFCEhIUbHEQIAX19f8vPz+fjjj5kwYQI5OTlGRxIWIMW9BrVs2ZKNGzdSr1497rnnHqPjCAFArVq1MJlMfPDBB0yePBk3NzejIwkLkOJeg8LDw1m1apWctQub06VLF9LT07l+/TouLi5GxxEWIMW9BqSnp1NYWMiyZcsA6Nixo8GJhPgjFxcXgoODmTt3LlprEhMTjY4kqqnc4q6UCldKJSml/q2Uml7GNmFKqX1KqcNKqXiLp7Rz7733HitWrODdd9/FZDKhVKWvahLC6vz8/EhJSeF///d/efrpp8nPzzc6kqiGO17nrpRyAY4B9wNpwB5ggtY68bZtPIGdwGCt9VmlVBOt9eVS9uW017lrrVmzZg0vvvgizz//vExaEjZr165dZGRksHv3bqOjiCLWWs/dHzihtT6ltc4HooGHS2wzEfhKa30WoLTC7uyUUrzzzjuYTCYp7MKm9e7dm0OHDpGQkGB0FFFN5VWaVkDqbc/PFn3tdj6Al1IqTin1q1JqkiUD2rMtW7bw7bff8t1333H27Fm6d+9udCQh7qh27dr4+/szY8YMTp8+zQsvvGB0JFFF5RX3ivRR3IBewFBgMDBTKeVT3WCOwNPTE09PT2bNmkVQUJBMnBd2wd/fn59//pmrV68SHh6Os7ZT7V151SYNaH3b89bcPHu/XSpwWWudDWQrpX4EegD/Lrmz2bNnFz8OCwsjLCys8ontiL+/Pzt27CAxMZGXXnrJ6DhCVEidOnXo1asXs2bNYuPGjUbHcTrx8fHEx8dXez/lfaDqys0PVAcB54Bf+PMHqvcCi7h51l4bSADGaa2PltiX03ygmpubC9z8FTc0NBRXV1f69+9vcCohKu7WIO2jR4/Spk0bzp49S5s2bYyO5ZSs8oGq1roA+AsQCxwF1mqtE5VSTyulni7aJgnYDBzkZmFfVrKwO5uNGzfy4osvcvDgQfbs2SODr4XdqV+/Pl27dmXGjBns2LGDGTNmGB1JVJIs+Wslubm5jBw5kuvXrzNo0CCj4whRaenp6XzyySckJyfTvHlzuT/DINa6FFJU0dmzZ4mLiyMgIMDoKEJUiZeXFz4+PrzxxhtS2O2QFHcLunTpElFRUQDMmDGD7t27y+BrYdeCg4NZvXo1GRkZxMfH8+mnnxodSVSQFHcLysrKwtXVlXPnzvH1118TFBRkdCQhqqV58+a0bt2auXPn0qJFCzp06GB0JFFB0nO3goiICPbt28fDD5e8mVcI+3P27FliYmI4f/48Hh4eRsdxOtJzN1hBQQEAV69eZe3atTL4WjgMb29vGjduzPz58wHIy8sjKyvL4FSiPFLcLeDy5cv4+flRWFjIm2++Sdu2bWnatKnRsYSwGJPJxKJFi8jPz+f1119n7dq1RkcS5ZC2jIWkp6dz11130bJlS8aNGyfzUYVD0VqzfPlypk2bxvPPPy/TmmqQtGUM5uXlxd/+9jeaNWsmhV04HKUUISEhLFiwQCY12Qkp7tUUExNDVlYW+fn5REVFYTKZjI4khFX4+vqSl5dXfDnksmXLSEtLMziVKIsU92owm83ExcWhtWbhwoXUrVtXBl8Lh1WrVi1CQkL429/+BoCbmxs5OTkGpxJlkZ67BZjNZlq3bk1YWBi+vr5GxxHCagoLC1m0aBHLli1j1KhRRsdxCtJzN9Ann3yC1hofH1nGXjg2FxcXTCYTc+bMKf5aYWGhgYlEWaS4V9GTTz7Jzp07MZvNMvhaOJUePXqQkpLC1q1bKSgooHv37qSnpxsdS5QgbZkqOnXqFM2aNeObb77hueeek8HXwqns3LmT69ev8/PPP3P58mWaNGlidCSHVdW2jBT3aurWrRvt27enV69eRkcRosbk5uaycOFCtm/fLiufWpn03GvI6dOn+f333wHYvHkzqampMvhaOJ1bg7Rfe+01ADIzM9m0aZPBqcTtpLhX0rZt24iJiQGQwdfCqfn7+7Nr1y6OHDlCTk4OGzdulGHaNkTaMlW0c+dOBg8ezEsvvYS7u7vRcYQwxJYtW2jSpAkbNmwwOorDkp57DQsLC8PFxUUGXwundu3aNaKiojh69Cht27Y1Oo5Dkp67lR0+fJjZs2cDcPDgQX755RcZfC2cXoMGDejatWtx7z02Nrb4sTCWnLlX0O+//86RI0cYMGAAw4YNIyMjg/vvv9/oWEIY7tYg7VOnTqGU4tq1a7Rv397oWA5D2jI15OTJk3Tt2pXnn39e5qMKUWTdunUEBQWxdOlSo6M4HGnLWNHVq1eLH0dGRsrgayFKMJlMfP7551y7dg24OSxe7lo1lhT3cly7do2AgAByc3O5cOGCDL4WohTNmzfH29ubuXPnAvDee+/xww8/GJzKuUlbpgLy8/Nxc3PjySefZO/evTL4WohSpKamsm7dOs6dOyeDtC1I2jJW5Obmxn/+8x+io6Nl8LUQZWjdujWNGjXivffeMzqKQIr7Hf3rX/8qnjQjg6+FKF9ISAiLFi2isLAQrTWRkZFkZGQYHcspSXG/g3PnzuHq6kp2djYrV66Us3YhytGuXTvc3d1ZtGgRSim6du2K2Ww2OpZTkp57BbzxxhusXr2axx57zOgoQti8xMREdu3axalTp2QZbAuQnrsF3f5DKD8/nyVLlhASEmJgIiHsR6dOncjNzWXVqlXFX8vLyzMwkXOS4l6KefPmFd+MsXDhQurUqSODr4WooFq1amEymYoHaaempuLv7y8rRtYwacuUIjMzk+zsbBo3bkybNm0IDQ2VwddCVMKtQdqffPIJI0eO5Nq1azRo0MDoWHZJ2jIWVK9ePZo2bcry5cspLCyUwddCVJKLiwvBwcHFg7SlsNc8Ke63ycnJYf/+/cDNvvu8efNk8LUQVeTn58fJkyfZtm0bACkpKezevdvgVM5Divttjh8/zuLFiwH48ssvycjIoEuXLganEsI+ubm5ERQUxMyZMwFITk7m0KFDBqdyHtJzL0O3bt1o164dvXv3NjqKEHbr1iDtuLg4/P39jY5jl6zWc1dKhSulkpRS/1ZKTb/Ddn2VUgVKqZGVDWFrYmNjOXPmDD169DA6ihB2rXbt2vTt21cGeBjgjsVdKeUCLALCgc7ABKXUfWVs9y6wGbC7BrXWmhdeeKF4aV8ZfC2E5QQEBLBz504SExMBeOWVV9i5c6fBqRxfeWfu/sAJrfUprXU+EA2UtiTiC0AMcMnC+WqE2WwmODgYT09Pfv75Z44cOSLtGCEspE6dOvTs2ZPIyEgAJk6cSLdu3QxO5fjKK+6tgNTbnp8t+loxpVQrbhb8qKIv2U9jvYiLiwsTJkxAKcWMGTPw9/fH3d3d6FhCOIzAwEC2bNnC6dOn8fPzk0sja0B5xb0ihfrvwKtFn5Yq7Kwtc+PGjeLHhw4dIiEhQT74EcLCGjRoQJcuXf7Qe7+14qqwjvKaymlA69uet+bm2fvtegPRRdeCNwGGKKXytdYbS+5s9uzZxY/DwsIICwurfGILe/HFFwkPD2f06NFERkbSu3dv7rrrLqNjCeFwgoOD+eSTT7h06RJ169YlPDycX375Rf69lRAfH098fHy193PHSyGVUq7AMWAQcA74BZigtU4sY/uVwNda63WlfM8mL4UsKCjAbDaTlpZGly5deO6556hfv77RsYRwSOvWrSM4OJioqCi01nKDYAVY5VJIrXUB8BcgFjgKrNVaJyqlnlZKPV21qLbF1dUVd3d3IiMj6datmxR2IawoODi4eJC2FHbrctqbmNLS0khKSmLQoEFcuHCBDh068N///d80atTI6GhCOLTo6GiGDx/O/PnzSUpK4sCBA4wbN87oWDZLFg6rpHPnzhXfCj1z5kx8fX2lsAtRA0wmE5988gk5OTm4uLiQm5trdCSH5LRn7rdkZGTQqlUrnnjiCZo1a2Z0HCGcwqpVq5gyZQqvv/660VFsnpy5V9Gbb77JPffcI4VdiBpkMpmKB2nfYosnf/bM6Yp7Tk4O//Vf/0V2djbZ2dmsWLECk8lkdCwhnEr79u1xc3NjyZIlAEybNo2YmBiDUzkWp2vLmM1mdu/eTXBwsAy+FsJAiYmJ/Pzzz6SkpHD+/HmaNWuGm5ub0bFsjrRlKqhWrVoEBweTn59PVFSUnLULYZBOnTqRk5PDZ599RqtWraSwW5hTFfeLFy9iNpsB+Oijj7jrrrsLC4QnAAAQlUlEQVRo27atsaGEcFK3Bmm/8847wM2e+2+//WZwKsfhVMV91qxZfPPNN5jNZj744AMZoSeEwbp27cqlS5dYv349+fn5vPrqq1y/ft3oWA7BqXrut95/+fLlzJw5k6efflqKuxAG27NnD2lpaezbt8/oKDZJeu4VoJRCKSWDr4WwIX5+fiQnJxMXF2d0FIfiFMU9OTmZL7/8Erg5+Prq1at07tzZ4FRCCLg5SDswMLD4hqb9+/cTFRVVzqtEeZyiuN+4cYOcnBwA5syZg8lkwsXFxeBUQohb+vTpw/79+/n1119p3Lgx3t7eRkeye07Vc9+yZQtjxozhxRdflPmoQtiYuLg43N3d2bJli9FRbIr03Ctg5syZMvhaCBvl7+/Pjh07SEpKAm7ecHj78gSichy6uOfm5hIYGEhmZia7d+/myJEj9OrVy+hYQohS1K1bl549e/Lqq68C8NRTT7Fx458GuokKcvi2zLFjx+jUqRMDBw4EIDQ0tMYzCCEq5tq1a0RFRXHs2DHq1auHp6en01/VJm2ZMnTq1InDhw+ze/du+vbta3QcIcQd3BqkPWPGDBo1auT0hb06HLa4nz59muzsbAAiIyPp1asXderUMTiVEKI8QUFBrFu3jsuXL1NQUMD27duNjmSXHLa4L1++nI0bN5KSksK2bdsIDAw0OpIQogIaN25Mx44dmTVrFgUFBSxZskSmNVWBw/fcJ06cSHJyMkOHDq3x9xZCVM2FCxf4/PPPOXfuHPXq1TM6jqGk516KixcvsmHDBoKCgoyOIoSohBYtWnD33Xczd+5co6PYLYcr7qmpqSxatAiQwddC2LOQkBCWLVtGbm4u8fHxxVObRMU4XHE3m814eXmRkZHBmjVrCA4ONjqSEKIKWrdujaenJwsWLKBt27Zyj0olOWzPfdq0aXz33XeMGzeuxt5TCGFZycnJbNmyhbS0NKddD0p67vzfeu3Z2dksX75cRugJYefat2+Pq6tr8SqRubm5xdPUxJ05THHXWhMUFERaWhrvvvsuTZo0kZXlhLBzSilCQkJ4//33MZvNjBgxgl9++cXoWHbBodoyqamptGjRglatWjF06FDatWtn1fcTQlif2WwmKiqKBQsWMHLkSOrWrWt0pBolbRlufgCzaNEiPDw8ZPC1EA7i1iDtt99+2+kKe3U4RHE/c+YMly9fLh58HRISImtSCOFAunXrVnzfSlZWFrGxsUZHsnkOsbB5bGwsubm51KlTh/z8fHx8fIyOJISwIBcXF4KDg3nzzTcJCQnhiy++4MEHH5STuDtwqJ67j48P3bt3p3v37lZ9HyFEzcvPz2fhwoV8/fXXhIWFGR2nxjh9zz0mJob09HS6dOlidBQhhBW4ubkRFBRUPEhb3JldF/esrCxee+01tNbMmTOH4OBgp73RQQhn0Lt3b/bt28dvv/3G559/ztq1a42OZLPsurjn5eXh4+PD1q1bOXXqFH5+fkZHEkJYkYeHB3379iUyMpIePXrIkgR34BA998DAQBo2bCjryAjhBLKysli0aBF79+7l3nvvNTqO1Tldz/3WLcgJCQkcPnyY3r17G5xICFET6tati5+fHzNmzABuFnvxZxUq7kqpcKVUklLq30qp6aV8/1Gl1AGl1EGl1E6llNUvVxkzZgy7du1ixowZ9O3bl9q1a1v7LYUQNiIwMJDNmzdz+vRpevfuzcWLF42OZHPKbcsopVyAY8D9QBqwB5igtU68bZsg4KjWOkMpFQ7M1loHltiPRdsyV65c4cyZM5hMJl544QWZjyqEk/n666/p3Lkz//znP/Hw8DA6jtVYsy3jD5zQWp/SWucD0cDDt2+gtf5Za51R9DQBsPqKXY0bN+aNN96QwddCOKng4GC++uorMjMzjY5ikypS3FsBqbc9P1v0tbJEAJuqE+pO0tPTOXHiBCkpKWzdulUGXwvhpBo3bkyHDh144403uHDhAj/++KPRkWxKRZYfqHAvRSk1AJgKlLqQ+uzZs4sfh4WFVekuswMHDrBp0ybS0tLo2rUr9evXr/Q+hBCOwWQy8dlnnzF+/Hh27NhB//79jY5UbfHx8cTHx1d7PxXpuQdys4ceXvQ8EjBrrd8tsV13YB0QrrU+Ucp+LNZzv3jxIu3bt+fJJ5/Ey8vLIvsUQtinNWvWMGLECObNm2d0FKuwZs/9V8BHKdVWKeUOjAM2lnjzNtws7I+VVtgtbdasWfj4+EhhF0JgMplYtmwZeXl5RkexKeUWd611AfAXIBY4CqzVWicqpZ5WSj1dtNksoBEQpZTap5Sy+KgUrTWvv/46p0+fZvXq1TJCTwgBQJs2bWjYsCHvv/8+s2fP5tdffzU6kk2wmztUCwoKWLZsGSdOnJDB10KIPzhx4gRbt27lyy+/pFOnTjRt2tToSBbj8Heourq6MnnyZFasWCFn7UKIP+jQoQMuLi4cOHDAoQp7ddhFcb+11MB7771H48aNZfC1EOIPlFKYTCbee+89zGazLEmAnRT3N998k8WLF7No0SI5axdClOree+/lxo0bLF68mL59+xafFDoruyjukZGR/Oc//8HDw4N27doZHUcIYYNuDdJesmQJe/fupVYtuyhvVmMX//Xu7u4sXboUk8kkMxOFEGXq1q0bFy5c4Pvvvzc6iuFsekC21prffvuNgwcPkpeXh6+vr9GRhBA27NYg7dmzZ9O+fXsApx29adNn7mlpabz11lvMmzePkJAQp/81SwhRvp49e3L8+HGio6NJSkoyOo5hbLpaent7M2nSJC5fvuy0P32FEJXj5uZGYGAgcXFxjBo1yug4hrHp4g4wZ84cTCaTDL4WQlRYnz592LdvH/v27TM6imFstuf+8ccf4+LiQkpKCsOHDzc6jhDCjnh4eNCnTx+mT59O8+bNiYqKol69ekbHqlE2e+bu5eVFVFQUQUFBuLm5GR1HCGFnAgIC2LFjB6Ghobi62ux5rNXYbHFv06YNSUlJMvhaCFEltwZpb9q0yaHH8JXF5hYOu7XNAw88QGFhYZUGegghBEBGRgZLly7l+PHjeHl52eVITodZOGzLli088sgj7Nq1C39/f6PjCCHsWMOGDencuTOPP/44zzzzjNFxapTNnbmbzWaGDBlCZmYmDz74YA0kE0I4sitXrrB8+XKSk5Np0aKF0XEqzWHO3M+cOcNPP/0kg6+FEBZxa5D2W2+9ZXSUGmVTxf3AgQNERkbStWtXGjRoYHQcIYSDuDVI+5tvviEjI8PoODXCZoq71pq//vWvbNiwgaCgIKPjCCEcSIsWLWjZsiXvvPMOqampRsepETZT3JVS+Pr64uvrK4OvhRAWZzKZOHbsmNMsQGgzxf3atWt8/vnnBAcHGx1FCOGA2rRpQ4MGDfjggw+MjlIjbKK4x8XFMXnyZLy9vWnevLnRcYQQDiokJIQPPviAv/71r0ZHsTqbKO6urq58//33MkJPCGFVHTp0wNXVlRs3bhgdxepsorjHxcXRrFkzWrdubXQUIYQDU0rRr18/YmNjHX7GquHFPT8/n8WLFxMSEmJ0FCGEE7g1SHv16tXk5eUZHcdqDC3uGRkZtG3bFjc3Nxl8LYSoEbcGab/00kv861//MjqO1Rha3OvXr4/Wmn79+sngayFEjenWrRuFhYU0bdrU6ChWY2hxX7VqFYWFhU5z3akQwjbcPkjbURlW3FNSUnj77bcxmUwy+FoIUeN69erF8ePHWbBggdFRrMKwqjp9+nTOnz9P165djYoghHBibm5u+Pv7M2/ePHJycoyOY3GGFffjx48TFhYmg6+FEIbx9/fnxo0bJCYmGh3F4gwp7lu3buXkyZP4+fkZ8fZCCAH83yDtyMhIo6NYXI0X98uXL/Pkk08SGBgog6+FEIYLCAggLi6O+fPnGx3Fomq8uCckJHD+/Hn69OlT028thBB/UrduXe677z42b95sdBSLqvHi/ve//52goCBq165d028thBClGjhwILt37yYtLc3oKBZTo8U9MTGRnTt3EhAQUJNvK4QQd9SwYUPuu+8+h+q9l1vclVLhSqkkpdS/lVLTy9jmH0XfP6CU6lnWvgYNGoSvry916tSpTmYhhLC44OBgPv/8c/bu3Wt0FIu4Y3FXSrkAi4BwoDMwQSl1X4lthgIdtdY+wFNAVFn7u3r1KmFhYdXNbIiUlBSjI1SLPee35+wg+Y1W0fxNmjTh3nvvZeXKlVZOVDPKO3P3B05orU9prfOBaODhEts8BHwKoLVOADyVUqVO3OjatSuNGjWqZmRjnDp1yugI1WLP+e05O0h+o1Umf79+/Vi1ahXfffed9QLVkPKKeyvg9mmyZ4u+Vt423qXtTEboCSFsWcuWLWnWrJlD9N5dy/m+ruB+Si7pWOrrWrZsWcHd2R5XV1e7vsLHnvPbc3aQ/EarbP7OnTvz7bffUlBQgKtreSXSdimty67fSqlAYLbWOrzoeSRg1lq/e9s2S4F4rXV00fMkIFRrfbHEvir6g0IIIcRttNaVXhO9vB9LvwI+Sqm2wDlgHDChxDYbgb8A0UU/DP5TsrBXNZwQQoiquWNx11oXKKX+AsQCLsByrXWiUurpou9/rLXepJQaqpQ6AWQBU6yeWgghxB3dsS0jhBDCPln8DlVL3vRkhPLyK6XClFIZSql9RX9eNyJnaZRSK5RSF5VSh+6wjU0e+/Ky2/JxB1BKtVZKxSmljiilDiulXixjO1s9/uXmt+W/A6WUh1IqQSm1Xyl1VCn1tzK2s9XjX27+Sh9/rbXF/nCzdXMCaAu4AfuB+0psMxTYVPQ4ANhtyQw1kD8M2Gh01jLy9wN6AofK+L4tH/vystvscS/K1wLwK3pcDzhmZ//fr0h+W/87qFP0v67AbiDEXo5/BfNX6vhb+szdojc9GaAi+eHPl37aBK31T8DVO2xis8e+AtnBRo87gNb6gtZ6f9HjTCARuLvEZrZ8/CuSH2z77+BG0UN3bp6opZfYxGaPP1QoP1Ti+Fu6uFv0picDVCS/BoKLfq3bpJTqXGPpqs+Wj3157Oa4F11d1hNIKPEtuzj+d8hv038HSqlaSqn9wEUgTmt9tMQmNn38K5C/Usff0lfoW/SmJwNUJMdeoLXW+oZSagiwHvC1biyLstVjXx67OO5KqXpADPBS0RnwnzYp8dymjn85+W3670BrbQb8lFINgVilVJjWOr7EZjZ7/CuQv1LH39Jn7mlA69uet+bmT8c7beNd9DVbUG5+rfX1W78+aa2/A9yUUl41F7FabPnY35E9HHellBvwFfAvrfX6Ujax6eNfXn57+DsA0FpnAN8CJScC2fTxv6Ws/JU9/pYu7sU3PSml3Ll509PGEttsBB6H4jtgS73pySDl5ldKNVdKqaLH/ty8nLS03pgtsuVjf0e2ftyLsi0Hjmqt/17GZjZ7/CuS35b/DpRSTZRSnkWP7wIeAPaV2MyWj3+5+St7/C3altF2ftNTRfIDo4FnlVIFwA1gvGGBS1BKrQFCgSZKqVTgDW5e9WPzx7687NjwcS9iAh4DDiqlbv2jnAG0Ads//lQgP7b9d9AS+FQpVYubJ62faa232UvtoQL5qeTxl5uYhBDCAdX4DFUhhBDWJ8VdCCEckBR3IYRwQFLchRDCAUlxF0IIByTFXQghHJAUdyGEcEBS3IUQwgH9f9NqCj+GhaK3AAAAAElFTkSuQmCC) 采用 [trapezoidal 方法](https://en.wikipedia.org/wiki/Trapezoidal_rule) 和 [simpson 方法](https://en.wikipedia.org/wiki/Simpson%27s_rule) 对这些采样点进行积分(函数积分为 2): In [35]: ```py result_s = trapz(y_s, x_s) result_s_s = simps(y_s, x_s) result = trapz(y, x) print "Trapezoidal Integration over 5 points : {:.3f}".format(result_s) print "Simpson Integration over 5 points : {:.3f}".format(result_s_s) print "Trapezoidal Integration over 100 points : {:.3f}".format(result) ``` ```py Trapezoidal Integration over 5 points : 1.896 Simpson Integration over 5 points : 2.005 Trapezoidal Integration over 100 points : 2.000 ``` ### 使用 ufunc 进行积分 `Numpy` 中有很多 `ufunc` 对象: In [36]: ```py type(np.add) ``` Out[36]: ```py numpy.ufunc ``` In [37]: ```py np.info(np.add.accumulate) ``` ```py accumulate(array, axis=0, dtype=None, out=None) Accumulate the result of applying the operator to all elements. For a one-dimensional array, accumulate produces results equivalent to:: r = np.empty(len(A)) t = op.identity # op = the ufunc being applied to A's elements for i in range(len(A)): t = op(t, A[i]) r[i] = t return r For example, add.accumulate() is equivalent to np.cumsum(). For a multi-dimensional array, accumulate is applied along only one axis (axis zero by default; see Examples below) so repeated use is necessary if one wants to accumulate over multiple axes. Parameters ---------- array : array_like The array to act on. axis : int, optional The axis along which to apply the accumulation; default is zero. dtype : data-type code, optional The data-type used to represent the intermediate results. Defaults to the data-type of the output array if such is provided, or the the data-type of the input array if no output array is provided. out : ndarray, optional A location into which the result is stored. If not provided a freshly-allocated array is returned. Returns ------- r : ndarray The accumulated values. If `out` was supplied, `r` is a reference to `out`. Examples -------- 1-D array examples: >>> np.add.accumulate([2, 3, 5]) array([ 2, 5, 10]) >>> np.multiply.accumulate([2, 3, 5]) array([ 2, 6, 30]) 2-D array examples: >>> I = np.eye(2) >>> I array([[ 1., 0.], [ 0., 1.]]) Accumulate along axis 0 (rows), down columns: >>> np.add.accumulate(I, 0) array([[ 1., 0.], [ 1., 1.]]) >>> np.add.accumulate(I) # no axis specified = axis zero array([[ 1., 0.], [ 1., 1.]]) Accumulate along axis 1 (columns), through rows: >>> np.add.accumulate(I, 1) array([[ 1., 1.], [ 0., 1.]]) ``` `np.add.accumulate` 相当于 `cumsum` : In [38]: ```py result_np = np.add.accumulate(y) * (x[1] - x[0]) - (x[1] - x[0]) / 2 ``` In [39]: ```py p = plt.plot(x, - np.cos(x) + np.cos(0), 'rx') p = plt.plot(x, result_np) ``` ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuczHX7x/HX5Uxko9JJ1F39UlTuTiLZKOfQ3VF0l3QTSSpRIcqhYpXSXblFUSJ0Ylch2ZwiwkaO1a7oQLLjnMP6/P74zmYas3bX2N2Znffz8ZjHfmfnMzOfx3fqcu31ub6fMeccIiISO4oU9ARERCR/KfCLiMQYBX4RkRijwC8iEmMU+EVEYowCv4hIjAkr8JtZZTObbWbfmdlKM+uaxbhXzGy9maWYWc1w3lNERMJTLMznHwAecc4tN7OywDdmNtM5tzpzgJk1Bc5zzp1vZlcDrwO1wnxfERE5RmFl/M6535xzy/3Hu4DVwBlBw1oAY/xjFgFxZlYpnPcVEZFjd9xq/GZWFagJLAp66ExgY8D9TcBZx+t9RUQkd45L4PeXeSYDD/sz/yOGBN3XPhEiIgUk3Bo/ZlYc+AB41zn3cYghPwOVA+6f5f9d8OvoHwMRkVxyzgUn1tkKt6vHgFHAKufcsCyGTQH+7R9fC/A55zaHGuici8pb3759C3wOmn/Bz0Pzz6dbYiIuPd07Tk/Hde6MS0vzfp+WhqteHZeSgmvf3rulpuI6d2bb3JUsoBZjErbQ68F07uJdrq25ixOL96B48UPEldjFBeV/5Zp/7qXxWd9y57kLaX/7Dh6sNovHqn9Gzwd89LxkGj0vmcbk1zZ775s5jwK6HatwM/46QFvgWzNb5v/dU8DZ/kA+wjk3zcyamtn3wG6gXZjvKSKxJikJ6tSBuDjvZ69e0KMHrFzp/WzeHMaNg8GDce+OY+OlzVny+hKWrCjB8nMWsKLSS/hGHuLCS5K54PVlnH/iFhoPqMvZo7rwUduTeOG57ZTs1d17rz59oP/LAcfvecdPXAj9P/COW18DDPTmMXCgN68oElbgd87NIwd/NTjnuoTzPiISg3IY7DOeH0JK7w+Zc+l/mXfDQOZdUwRXMZUrB63kypPTeODlS6nx3xZU+WoERU4qD93f9F7/vL0wux+zb72VksnTISHB+/2wYdkfz58PzZp5QT/zOJoU9J+JAX+yuGg1e/bsgp5CWDT/gqX5B0hMdC493TtOT3euc2fn0tK836elOVe9unMpKe6nu59yb/Te6G5loqtQ/oC7MO5X1/HWre6ds55wP87d5A5tS3eufXvvNmGC99zOnb3j9HTvlpjozX/q1L+Oo40/buY63poLo050PJmZi5S5iEg+Cszsfb6/Z/bVq0Pz5rh3x7Gk/6d8FNeOxFG/8UvcRTQ+ZSkN7zmNBu+048zpow+PHzwYrrsOGjXyXj8zI/f5ojM7Pwozwx3D4q4Cv4gUrMxgn1kr37Dhr2C/8NkZTCh1Lx++t5cTqp7CzSd9SYtH/sFVL9xC0aQpfw/2mc8vhAE+Kwr8IhI9ArN88IJ19+5QuzarZ25ibJkHGD96D6XPOY3WJ3/ObU+dT7U+t0JiYswH+0DHGvi1O6eI5I+kJC9Aw+HF2g0bICmJHTvgjRV1qNX+IhrMepKMpSl8MnE/q8pcydOTalCt+Pde0B882HtulSqHF1bBC/4xFvTDoYxfRPJONvX7JQ168MYlr/NBUkkanLGa9o9X5MbXbqZY0ifK7HNApR4RiTwh6vf7m7ZiUttPGP5yBr8VPZOOv/Sl3S07Oe3NAV5QV7DPMQV+ESl4wbV78Mo5Dz7ItgGv8cYDy3n1xyZc9PuXPDTsPJqvSaDolf+EBQu8HvnAmr+CfbZU4xeRgnGU2j0+Hz/1GUnXE9/mvJplWV/pWqbXH8znKafS8s2bKPrE43DffV7Q79Xr8OuoZp+nFPhFJDyZwd7n8wK2/6radRn/4L5rVlNz6jOUWr2MlTN/460f61HjhbawcePhxdrM5wUu1kqeUqlHRHLvKO2Y62Zt5FlfV6ZPO0iXThk8tG8oFYb2Uv0+D6jUIyL5JzDL9/tx64nc074odT7txYVb5/HDt3vou6E9FZ7ucrh0ozbMiKCMX0RyJoss/7eLGzBgZCXG/3gVD7XbzSNftqL8pxO8IB/c1SPHlTJ+ETn+jrJwu2sX9P26GRc/2pASbh9rVmbQr/kSL+irdh/RlPGLSNZC9OFnNGvB6FuS6JtQlvqnr2Fghw1UWTdT7ZgFQH38InJ8HGXhdvaUHXRb3ZG4dV8ztMksrnjv0b9flauSTr5SqUdEjo8QC7dpf5TjlvbluW9xZ/qcM47kN3/gijN+OfwclXSiijJ+Eckyy//zyrok/O9EXvruRrq130X35OaUnvaBFm4jhDJ+ETl2IbL8GannU+OB2nzze2W++foQfZp+4wV9LdxGPWX8IrHqKO2Zj7x6Lot+P5dX71tG010TtXAboZTxi0juBGX5hw7ByJXXUOPRG6hSegsr1xSn6bCG2kenEFLGLxJLssjy1/+jMR2GV2fPwRKMfDCFSzYmKcuPAsr4RSR7QVl+RgYkLGvANU/F0/KMJSxYXYFL+t6sLL+QU8YvUthlkeWvOacJ7V6+jFLFDjDq4RWcu366svwoo4xfREILkeUPXVafa3vX4+6qc5m18jTO7XmbsvwYooxfJBb4e+5T73ySe27bA8Dbj6Qoy49yyvhF5LDAzdUAVz6O0ac+wVXXlaTlmUuY/d2pyvJjmAK/SGEUUN7ZuhVuabKHl5/fw+wnZ/JYzS8oWtQ/ThdhxSSVekQKixCLuDNbj6bdwg60LjaJAfPrU/ICbbVQmKjUIxLrArL8/fuhx9OlaDfnXsb4WjJk1uVe0Adl+aLALxLVAmv5/oD+fYfB1Pm/31nzyRqW3/wsDVJHwYgRf6v5q5Yf2xT4RaJZUKvmhA+KU3vKE/w77Vk+qTuUk1/tB1Wrehl+0CZsErtU4xeJNiFq+Xu7PUm3tIf5YtlJTLz+DWq2qAwLFqhVs5DTN3CJxIqgxdl1S3Zw240+LvbNY8Qdsyn3xhB9K1aM0OKuSKzIXJzt1YuJr27h2rrQ+eJkxr35J+XKuiPHaRFXgijjF4kGQeWd/fuh+71bSRq/nUlN3uKf73VXlh+DCizjN7PRZrbZzFZk8Xi8mW03s2X+W+9w31Mk5gQs4m7aBPHX/Ela4gq+eWEW/zzjt8PjlOVLDoSd8ZtZXWAXMNY5VyPE4/HAo865Ftm8jjJ+kUAhFnGT24zkrnmdeaj0SHp+dTNFztEFWbGswDJ+59xcID2bYbmemEjMC8jynYMX/1uSO7/4D2N3tOTJGfW9oA/K8iXX8mNx1wG1zSzFzKaZ2UX58J4i0c8f0Hf3eIY2rXbx7ktbWNRiEDekvqkLsiQsxfLhPZYClZ1ze8ysCfAxcEGogf369fvrOD4+nvj4+HyYnkgECSrvpPniaDVnEJeuncj8O7+k9Osv/q2rR+Wd2JKcnExycnLYr3NcunrMrCowNVSNP8TYVOBy59y2oN+rxi8SUK//Ymkcd92RwZPFhtD1qbJYynJdkCV/U6AXcB0t8JtZJWCLc86Z2VXAROdc1RDjFPglNgVl+S7dx6tNpzHwu5a8V7Er9ZOfhipaxJUjFVjgN7PxQD3gZGAz0BcoDuCcG2FmDwKdgIPAHrwOn4UhXkeBX2JTQEDfXyaOzvfvY9HULUzxXcc5KZ/AJZf8fayyfPHTlg0i0cznY8sjz3HLymeouGU179R7k3LPPg5DhijDlyxpywaRaBL01Yjf/hTH1TP6U29JAh/WfYlyrwzUrpqSZxT4RQpCQI/+1KnQ4PoMBhbpw4DhcRQpVeLwOPXoSx5QqUckv4RYxB3aeAYvrW3KhxU7cPUXz2kRV3JFNX6RSBcQ0A+c4C3ifj1lM1N913J2SqIWcSXXVOMXiXT+sk36YwNocv1efp29hnk3vcDZqXN0Ja7kKwV+kbwUtIj747Y4as8eQPX5I/ik7lAt4kqBUOAXyUsBi7gLF0KdWgfp8udQhg0vRtHSWsSVgqEav0he8/mYfOdkOn19L2+f2JVmX/bUIq4cF1rcFYkUAd07zsHQoTBsyH6mbrmamiljtIgrx40Wd0Uihb+8c3Crjy5dYMyb+1kQ18wL+lrElQigjF/keAjq0d/9s4/WtdPYW7oCk4vcQflPJ6i8I8edMn6RghSwiLtlC9RvWZaTiu0kae15lJ8wwgv6oEVciQgK/CLHgz+gr39wGLWvOkCjjE95O34MJVLXqbwjEUelHpFjFVTeWbQIWjU/wDNbH6RD270wfLj3mMo7kkdU6hHJbwHlncREaN40g5FlutFh+CVQsuThcSrvSIRRxi8SDp+PN29Oos/K2/n4xHu00ZrkK/Xxi+SHoB79/v3h7f/t47Ofa3BBymT16Eu+UqlHJD/4yzsZf/jo3Bk+nriPBWUbeUFfi7gSJZTxi+TSn7/5uKt2GjtO/Qcf7riREz99X+UdKRAq9YjklYDyjs8HLVvCGWV8jPnsVEqkLFF5RwqMSj0iecVf3vll9Xauuw4uO3cH4zZe5wV9lXckCinwi4QSuI9+XBzr2j1Hnav2c1fVBQxbXIciSVO9TF/76EsUUuAXCSWgR/+bbyD+prL0rv4JT0ytg703TlswSFRT4BcJxR/QZ90zliaNMnjt4tdof/FCSE1VeUeingK/SKagr0mc/Hkcred2ZtIf9Wl1+iJISNDXJEqhoMAvkimgvDNiBHR98CAzKtxJveG3aQsGKVTUzikSwKX7GNR4DqN+bsSMMjdz3szX1aMvEUt9/CK5FbS75qFD8FiXP/l8yh5m/Hwxp6dMV4++RDT18YvkVkBp58ABaNdmP19/sIk51/fj9NSvtIgrhZYyfoltPh97e/bjjh8GcXDtD0yq/zonvDxI++hLVFCpRyQngso727dDiwa7OOObqYwZsY8St7f6e5BXeUcimEo9IjkR9N2419fZR/XURMYtu5gSKYuPHK/yjhRCyvgl9vh8/NQ1gYZze3Pb7jE8+3VjrKo6dyT6qNQjkpWg8s7atdDw+v10+7Unj6S0U+eORC2VekSyElDeWboU4usepF+R/l7QV+eOxKCwAr+ZjTazzWa24ihjXjGz9WaWYmY1w3k/kRwL2l2TgQOZ02YEjevt4bXSj9Fu/v3aXVNiVrgZ/1tA46weNLOmwHnOufOBDsDrYb6fSM4EZPkASdOLceusTozfdRM3T22v3TUlpoUV+J1zc4H0owxpAYzxj10ExJlZpXDeUyRHMgN6r16Mf+V32rfLYGrD4TRIHaXyjsS8vK7xnwlsDLi/CTgrj99TYlXQ7prExfH6SU/y+MP7+LzB81w99kHtrilC/izuBq84q3VH8kZAecc5GNQjnYTBjjl9v6B6pd8Pj1N5R2JcsTx+/Z+BygH3z/L/LqR+/fr9dRwfH098fHxezUsKI39Ad0/1oseefnw2wce8+aU5/cp/g6/F33v0Vd6RKJScnExycnLYrxN2H7+ZVQWmOudqhHisKdDFOdfUzGoBw5xztbJ4HfXxS+4F9ehnZEDH27ax4qP1fDqnLBXqXnx4rHr0pZApkD5+MxsPLAD+z8w2mtl9ZtbRzDoCOOemAT+a2ffACKBzOO8ncoSA8s6+fXBni92kzVjHrK9OoMKE17SIKxKCrtyV6OfzsbvHM/xr5bOUWfk14xefR6n/0xYMUvjpyl2JHUHdO+kujoaLB3DaVx8yKfkUL+iDFnFFsqDAL9EnoLyzeTPE197PlRs/5K1lNSk2Sj36ItlRqUeik89H2kNDuXFOb+7+8036LGqOdtiUWKPdOaVwC+reWbUKGtXfz+Obu9M15X7tsCkxSTV+KdwCyjuLF0P9ehkMKtbXC/ragkEkV5TxS/Tw+Zh979vcMacLb5Z7hBZzunubram8IzFKpR4pfILKOx9/DB3uO8jE9BuIT3lF5R2JeSr1SOETUN55+23o1OEgn1Zo4wV9lXdEjpkyfoksQVk+Ph8vNprOy+ubMP3ktlw4c7jKOyJ+yvilcAjaYfOpZ0vxvzV1mZtenQsnD9AXqIgcB8r4JfL4fGQ82ZvO259j6cytTGv4MqcM7AZDhijDFwmgxV2JXkHlnX374O4WPrbOWMond4yn3BtDvMdU3hH5G5V6JHoFlHd27YLmN+wlY+Fipr24lnJlA5IBlXdEjgtl/BIZfD62PjqIZov7ccnP03hjyZUUPVeLuCJHo1KPRJeg8s5PP0Gj+D9pmTqM55Y3xS5Vj75IdlTqkegSUN5ZtQquveYgHfa+wvMpTbH/qUdfJC8p45eC4/OxsP1IWn3ZjYTST9N23gPq0RfJBZV6JPIFlXemTYN72x7k7fQWNE15XlswiOSSSj0S+QLKO2PHwn33ZDDlpHu8oK8tGETyjTJ+yVtBWb5L95HQaCavft+Izyq2pdrn2oJB5Fgp45fIFJDlHzoEj/YqxZh11zA//SKqfaAtGEQKggK/5C1/QN/3RF/atNrFNx/+xNybBnNW6jyVd0QKiAK/HH9JSX8L6NstjqbLBrBv6gymX/88Jw1/FqpW9TJ8/18DIpJ/FPjl+Aso7/zyC1xXaz8Xfp/EpJd/pfQJAf/JqbwjUiC0uCt5w+djdadXaDLnCTpmvMYTC2/GqmoRV+R4Uh+/FKyg7p25c+HWVgcYvO1+7kl5TD36InlAXT1SsALKO5MmwS03Z/BOuQe9oK9FXJGIooxfjl2IHv2XGk/nxXXNSKx4L5fNGqoefZE8pFKP5L+AgJ5RLo5unfYxe+LvJG2vQ5WUqSrviOQxlXok//m7cvb06MctTfawKvFH5t30AlVSv1R5RySCKfBL7gT16P/2ZxzxCwZRfuYkPo1/gbjh/dWjLxLhFPgldwIWcb/7DmpdcZBmv7/F26/spESZYofHqUdfJGKpxi+55/Px+b/Hcte8zrx4Qh/toy9SQFTjl7wTVN7538Q42szvxKT0BrRNaq2N1kSijAK/ZM9f3sn4w0f37pDw/AHmVWxFvZThWsQViUIq9UhoQT36uzb5aFMnlR0lT+GDondQ4bP3VN4RKWAFVuoxs8ZmtsbM1ptZzxCPx5vZdjNb5r/1Dvc9JR8ELOL+9BPUaVyOU0v4mL7+XCq8/7rKOyJRLKzAb2ZFgVeBxsBFQGszqxZi6JfOuZr+24Bw3lPyiT+gL2w/klpXHuTespP533XjKJG6TuUdkSgXbsZ/FfC9cy7NOXcAmAC0DDEu13+KSAEIWsQd80kcLWY/wsgtLXnk/ERsaIJ69EUKgXAD/5nAxoD7m/y/C+SA2maWYmbTzOyiMN9T8krQIm7/vgdIrngLzYY3gZIlD49TeUckqhXLfshR5WQ1dilQ2Tm3x8yaAB8DF4T5vnK8BC7ixsXh6zGI1pf9xP4yO1lUujUVPxsXehFX5R2RqBVu4P8ZqBxwvzJe1v8X59zOgONPzew1M6vgnNsW/GL9+vX76zg+Pp74+PgwpyfZylzEHTiQVb/E0apFWZqc8AsJa/9J8ZRvQi/iKuCLFIjk5GSSk5PDfp2w2jnNrBiwFmgA/AJ8DbR2zq0OGFMJ2OKcc2Z2FTDROVc1xGupnTO/BLVq4vMx5bZ3uH/R/bxw2QTaXTAfeveGIUPUpikSwQqkndM5dxDoAkwHVgHvO+dWm1lHM+voH3YrsMLMlgPDgDvDeU85DgJaNQ8dgr7PleLBBXcxdWc87ap8AQlaxBUpzHQBV6zy+fB1H0Db1GfZ8d0mJtZ/g9NuqA4LFniBP+CvAZV3RCKTvohFji6ovPPtt3DLTfto8tMIht61lOL/HeY9pitxRaKGNmmTowso77zzDjS4PoN+h/ryyvAiFC+t7ZRFYoky/sIsKMvft9nHI/HL+PyPmnxQ/j5qfP6S9tsRiWLK+OVIAVl+airUaVyWzbvLsvj3KtT4oJ/22xGJUQr8hU3gtgv+gD7ltne4+tK9tC05mck3jqB8aor22xGJYQr8hU1Alr9/PzzaqxQPfdWaT3bWp9v5SdpvR0QU+AuFEFl+asfnufbC3/khcTXLbu7PNaP+o/12RATQ4m7hELQ4+/7o3TzUOYNe+/rQtc027NXhatUUKYTUxx9rQmy7sPvhp+j6Yzfmfnsi468fyeUtztQFWSKFmAJ/rAnK3pfM3kmblruovfMzht8xn7JvJCjLFynk1M4Za/w1+ownezPo8XSaNnU8c/kU3hrlKFv2yHGq5YtIJmX80SK4tAP8uHwH99y+h2LrVzP25o+pPPoZZfkiMUQZf2EU2K2T2aa5YQMuMYmRL+/hqlpGy5PmMGtkKpUr7D78PGX5InIUyvgjWXDmvmEDvzRqR4dTP+bXVemMrTeKi0c9qixfJEZpcbewCNGtQ/fuuGtq8857Rem+rA2d0gfSa0QVStze6u9BXh07IjFFpZ7CIuDK20ybtpXhpvtP5cWfb2fGjUN4JvUeSqQsPvK52nZBRHJAGX8kyCLLP1SrNiPeKsHTS1vyULvdPPFlE0pM+1g7aooIoIw/uoXI8ldvqUj8f87jnXVX8+Vsx9PNvvGC/uDB3jgt4IrIMVLGX1CyyPL3XlGXgW9UZMTaevT910o6lRxN0ReH6MpbETmCMv5oEyLL/yz1Amp0qsO69FNIWXqILuOu8YJ+4DjV8UUkTMr481MWWX7aBQ155LXzWLHtLF65bzlNd03U/joiki1l/JEqi4uwSEpizx54ZnFTruhZn8vjfmTluhI0HdbQC/rK8kUkjyjw57XAkk5cHPTogWvWnPe/q061CzJYtbMy37wwi95XfEapUv7naOFWRPKQSj15IYuSDrVrsyBxG91/7MSfKWt5uWESdd9/SFfeisgxUamnoB2lpAOw7veTuLX9idz5VVceqDiZJSOXU7fyhsPPV5YvIvlEGf/xEmJfHZo3Z9OLE3n2oS18tOGfPNphF90+v4nS0z7QRVgiEjZl/AUhxHfd0r07jB7Nb/3e4NHLvuDShqdSodgO1q46xJMNl3pBXxdhiUgBUuDPrWxKOpt9JenefhsXfdifjJSVrBw6g+drfUKFk5zXmVOlyt+DvTp2RCSfqdSTE4GLtZnlmR49YOVKqF4dmjfnp4SJDH7kV9774WruunkvTy6/gzOnj1ZJR0TyjEo9eSlESybNm0Plyqzo+S73nvslNRufygm2h1XfHuTVuxd5QV8lHRGJQMr4s3KUlsxDS5Yys8ajvNR5Hd+e0oCHzp1Gx7t2UuHbZF1xKyL5Rl/EcjwcraRTpw47Oj7O2IklebVqAqV2buXhro673m9JyWkfqaQjIvlOpZ5jldVi7fz5f5V0luy4gP9c/S1VPh7GlzcOYKR1ZNmSDNpd/q0X9FXSEZEoEpuBPwfBfkuZqgxrs5iaB77mtruKcU6Z31iVcpBJD8+n7uxnsSGDveeqS0dEokzhLvUElm4yjwGmT4c5c47ozNk1cjxT+ixm/P5bmDvnEC2aO+4tMpb4m8pRZOEC1e9FJKKoxp8pqzr9woUwc6Y3JiEBtm+H5s3ZMWI8Sf2+5iP7F9NnGHWuNVqX/IhW/S+nXIfWkJio+r2IRKQCq/GbWWMzW2Nm682sZxZjXvE/nmJmNY/5zQJLNJnHPh/063f4eNeukKUbqlX762V+SNnF8LaLaHziAs6qczbv/nkrjfYn8kPyJqb56nD3qHjKbdvgBX3V70WksHHOHfMNKAp8D1QFigPLgWpBY5oC0/zHVwMLs3gt59LTvVvfvqGPJ0xwrnNn59LSvOP27b1bWtrh4/R073716s6lpDjXubPb+uVKN5l/uQfabHfnV93vKvGra9fyDzep0Ui3ff4Kb2xamnOJid7Pzp2913HO+5mY6EREIo0XwnMfu4uF+e/GVcD3zrk0ADObALQEVgeMaQGM8f8js8jM4sysknNu8xGv1r2797NPn9DHCQlQq5aXwY8bd7h0E1AiyvjDx5re77L49mksuPRT5p07lE0jHdfGv0OD9Yk8cNkqanz0L4rc3cbL6FeuPJzZZ5ZxMjP7Zs20WCsihU64gf9MYGPA/U14WX12Y84Cjgz8mQJr/YHHPh8MGQLjxuEuvZQtX29gXWpxVp7zLCvavsiKNSVYft4BTjv7ca5YsoLaPZvRafLd1FidQLGK5aH7DO911q5WsBeRmBVu4M/pamzw4kPI5/UrW5ZDh2DfOdW47OWxVDuvDlvOuZfNQ8exeWtRNp0zlp/qD2LDDdv5vlwGxevu5vwTt3Dx7YOpsT6RW0/5mZpzm3BSpzv92fxy6JTgBfjrrvP+YgAvwAe2YSrYi0gUSE5OJjk5OezXCaurx8xqAf2cc439958EDjnnXggY8waQ7Jyb4L+/BqgXXOoxM1fUDuIwTiznKH8onfIl9nJKtZOp9FsKlUrt4KzG1Tl7YgJnD3+cf9QoQ8VBj3lP7tMH+vf3jm+80SsHBWbzar0UkUKoQNo5zawYsBZoAPwCfA20ds6tDhjTFOjinGvq/4dimHOuVojXcvvbdaCYZWBPBwTyrIL6dddBo0be74cNg27dvOPMAK9gLyKFXIH18ZtZE2AYXofPKOfcc2bWEcA5N8I/5lWgMbAbaOecWxridZxLT/fuBAZyBXURkZB0AZeISIzRJm0iIpIjCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYUO9YnmlkF4H2gCpAG3O6c84UYlwbsADKAA865q471PUVEJHzhZPxPADOdcxcAs/z3Q3FAvHOuZmEN+snJyQU9hbBo/gVL8y9Y0T7/YxFO4G8BjPEfjwFaHWWshfE+ES/a/8PR/AuW5l+won3+xyKcwF/JObfZf7wZqJTFOAd8bmZLzOw/YbyfiIgcB0et8ZvZTOCyal8SAAAEA0lEQVS0EA/1CrzjnHNm5rJ4mTrOuV/N7BRgppmtcc7NPbbpiohIuMy5rOJ1Nk80W4NXu//NzE4HZjvnLszmOX2BXc65oSEeO7aJiIjEMOdcrkvpx9zVA0wB7gFe8P/8OHiAmZUBijrndprZCUBD4JlQL3YskxcRkdwLJ+OvAEwEziagndPMzgBGOueamdm5wIf+pxQDxjnnngt/2iIicqyOOfCLiEh0ytcrd82ssZmtMbP1ZtYzizGv+B9PMbOa+Tm/7GQ3fzOLN7PtZrbMf+tdEPMMxcxGm9lmM1txlDGRfO6POv9IPvcAZlbZzGab2XdmttLMumYxLiI/g5zMP1I/AzMrZWaLzGy5ma0ys5BVhwg+99nOP9fn3jmXLzegKPA9UBUoDiwHqgWNaQpM8x9fDSzMr/kdp/nHA1MKeq5ZzL8uUBNYkcXjEXvuczj/iD33/vmdBlzmPy4LrI2y//5zMv+I/QyAMv6fxYCFwLXRcu5zOP9cnfv8zPivAr53zqU55w4AE4CWQWP+uijMObcIiDOzrK4PyG85mT9E6MVqzmuhTT/KkEg+9zmZP0TouQdwzv3mnFvuP94FrAbOCBoWsZ9BDucPEfoZOOf2+A9L4CVx24KGROy5hxzNH3Jx7vMz8J8JbAy4v8n/u+zGnJXH88qpnMzfAbX9fypOM7OL8m124Yvkc58TUXPuzawq3l8vi4IeiorP4Cjzj9jPwMyKmNlyvItNZzvnVgUNiehzn4P55+rch9POmVs5XUUO/lcrUlafczKPpUBl59weM2uC1+J6Qd5O67iK1HOfE1Fx7s2sLDAZeNifOR8xJOh+RH0G2cw/Yj8D59wh4DIzKw9MN7N451xy0LCIPfc5mH+uzn1+Zvw/A5UD7lfG+1f1aGPO8v8uEmQ7f+fczsw/yZxznwLF/W2v0SCSz322ouHcm1lx4APgXefcEde9EOGfQXbzj4bPwDm3HUgCrgh6KKLPfaas5p/bc5+fgX8JcL6ZVTWzEsAdeBeBBZoC/BvAzGoBPnd4P6CClu38zaySmZn/+Cq8dtlQtbhIFMnnPluRfu79cxsFrHLODctiWMR+BjmZf6R+BmZ2spnF+Y9LAzcCy4KGRfK5z3b+uT33+Vbqcc4dNLMuwHS8xYlRzrnVZtbR//gI59w0M2tqZt8Du4F2+TW/7ORk/sCtQCczOwjsAe4ssAkHMbPxQD3gZDPbCPTF606K+HMP2c+fCD73fnWAtsC3Zpb5P+1TeBdARsNnkO38idzP4HRgjJkVwUt233HOzYqW2EMO5k8uz70u4BIRiTH66kURkRijwC8iEmMU+EVEYowCv4hIjFHgFxGJMQr8IiIxRoFfRCTGKPCLiMSY/weoVZxsAST89wAAAABJRU5ErkJggg==) ### 速度比较 计算积分:$$\int_0^x sin \theta d\theta$$ In [40]: ```py import sympy from sympy.abc import x, theta sympy_x = x ``` In [41]: ```py x = np.linspace(0, 20 * np.pi, 1e+4) y = np.sin(x) sympy_y = vectorize(lambda x: sympy.integrate(sympy.sin(theta), (theta, 0, x))) ``` `numpy` 方法: In [42]: ```py %timeit np.add.accumulate(y) * (x[1] - x[0]) y0 = np.add.accumulate(y) * (x[1] - x[0]) print y0[-1] ``` ```py The slowest run took 4.32 times longer than the fastest. This could mean that an intermediate result is being cached 10000 loops, best of 3: 56.2 µs per loop -2.34138044756e-17 ``` `quad` 方法: In [43]: ```py %timeit quad(np.sin, 0, 20 * np.pi) y2 = quad(np.sin, 0, 20 * np.pi, full_output=True) print "result = ", y2[0] print "number of evaluations", y2[-1]['neval'] ``` ```py 10000 loops, best of 3: 40.5 µs per loop result = 3.43781337153e-15 number of evaluations 21 ``` `trapz` 方法: In [44]: ```py %timeit trapz(y, x) y1 = trapz(y, x) print y1 ``` ```py 10000 loops, best of 3: 105 µs per loop -4.4408920985e-16 ``` `simps` 方法: In [45]: ```py %timeit simps(y, x) y3 = simps(y, x) print y3 ``` ```py 1000 loops, best of 3: 801 µs per loop 3.28428554968e-16 ``` `sympy` 积分方法: In [46]: ```py %timeit sympy_y(20 * np.pi) y4 = sympy_y(20 * np.pi) print y4 ``` ```py 100 loops, best of 3: 6.86 ms per loop 0 ```