{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 符号积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "积分与求导的关系:\n", "\n", "$$\\frac{d}{dx} F(x) = f(x)\n", "\\Rightarrow F(x) = \\int f(x) dx$$\n", "\n", "符号运算可以用 `sympy` 模块完成。\n", "\n", "先导入 `init_printing` 模块方便其显示:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from sympy import init_printing\n", "init_printing()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sympy import symbols, integrate\n", "import sympy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "产生 x 和 y 两个符号变量,并进行运算:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAAAlCAYAAADY4B6YAAAABHNCSVQICAgIfAhkiAAAA4FJREFU\naIHt2UuIHFUUxvFfEskYoxFfmBiNY4JEGYzGiEqiaNABSVw4IrNQcZOVulBIFAV3IkRRwfiAgIsy\nii9ERfCBJMYEUUQFdaGoEB8Iig9E1BAFExenm+kuu0xN1e2Zsak/NN236ta5H7fOveec2zQ0/J+Y\nlcDGgQQ2Gnrw7HQLmInMrvn8cnyeQsigUXdi1+PlFEIaunlG/ZczkNSZlCPwB/Yn0jJQ1JnYUWxP\nJaRhgq04ZrpFzFTqeOzR+DmVkEGj6sSuxIcphQwah1R8LlWadR5W40iswZ3YncBuSqZU43MJbByO\nzR3tcezF4gS2UzGlGo8VgasuK0SqtqzVXiDOHcYT2E7FlGq8FlcmsDNLLLP2QdCIEL0yge1UTKnG\nTBQHqXkc9/XBbkpKa+yVFSzBnIL+c3AYfqumq5AN+A6batpZoXpAPhi1NB6K94X792INNlbTVcjl\nQnR7/OEatrKazxcxaY15j90oIt4lBf1Tn2ZdhONbNhfiMixKaD8FtTWehLOwDW8U9HmhpK1VeAD3\n43lRpd2Ou/GEiLJLxZZyIPdZMBnROTLlPXYYD4oJuzp370ZxDpJU43XYh3m56ydiS4nnl+JhE6sh\nE4fhq8VWsl/67aRNpvzEPoK5uAkf5+69i6frCOkVvHZgSExCJ+vwagmbm3CbiePE+fgFb+Mb4cVZ\nBa0puQBv4S+xtDv/BZmPs7GrHwN/qrvigCfFxn0wTsm1v8VdKUSVIFPOYxcK51mMvzHWcW9ULPeR\nOkKKUpPtugPYXOHd+0rY/LLj93IhfmcldcU8hjN7XF+Cc4Un5tmAD1q/v299j+N3vNLR70L8hE+S\nKM1xhXiTR7Xao7ihgp3r8afIfdssK+ibgszk0q3X8GLu2psi4Nai6Nhwp1gOa1vt9brfahHzcA/O\naLVHRWDY2zHeLZWU9oeT8VlHe0jk8LVPr4om9ldRKFzaag/jqxL21omJG8FpOFV4bJs7xDKeKXwt\nUsE2m0UcqR24/qv824GrxOR8UdLeLrEcV4nIer5Ia7aKfe8lvFNRaz+4GY/iIbGqzhFO9VE/B10r\ntoN7cXE/B0pIpnpJO1sEtW2pxBQxJN7iD/p3sJGaLSKVKsNTuguDMbGqTk8tqhevS/NvwUzkR1HS\nwgnYg2tSGS86HmyzCO/5d8k3COzBcaLyGsOtylWWDQ0NDQ3TzD+oDaLToDx6/gAAAABJRU5ErkJg\ngg==\n", "text/latex": [ "$$\\sqrt{x^{2} + y^{2}}$$" ], "text/plain": [ " _________\n", " ╱ 2 2 \n", "╲╱ x + y " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y = symbols('x y')\n", "sympy.sqrt(x ** 2 + y ** 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于生成的符号变量 `z`,我们将其中的 `x` 利用 `subs` 方法替换为 `3`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEoAAAAlCAYAAADlcn/+AAAABHNCSVQICAgIfAhkiAAAA49JREFU\naIHt2WuoVFUUwPFfXvEmlqI9vGbqTYmIwDKiKIuSuiDZh4y4HzJKECIqSJDyY0RBEtkHiyCImAx6\nEdWH3mgvougFWVBkIAUFYS+kMrHMPqwzdGY65845c851Jps/XIZ99tlrr1mz9tprrcuAAXVyRA0y\nDtYg43/BU71W4FAwpeL6U7CzDkX6naqGWoUX6lDkcOdJ1Y39n6DKlzwav+GvmnTpa6oYagzb6lLk\ncOZBHNNrJQ4VVTxqDn6sS5F+p1tDLcPHdSrS70ztcl0dacE5OA+zsBx34K2KMvuOpyuuPwqbUuNx\n7MX8inL7imNFIK/CUpFWLEnGM0XNOF5R7qTRzdFbiVcq7vspzsWuZLwg+fyyotyiLMZt2IcDmIGN\n+K7OTRoi2ayTR7G5Zpl5nIQfcEXq2dXYoaTjLMRQztyQ+rsF63C36i2fpYp90efwk9Ybf7rwrrVF\nNzsSH4obKYvl2FBUWAEuE4Zq7j1aQVajwPpp+APvZ8x9gefzFrbnURvEzXNxzvt1dgsuxNxE3oiI\nffNqkp3HHOF1+zLm9uDsIkIW4AxsxWs57zxbQM4o7hMGuKpt7kZRHy7GL+KmS//NLKJoDg2dPWpI\nFPJZHvVNokPhOHWNsPj0tucnYkuB9Q8IF78Zn7TNvYcniipSkoZiR/dhEczTMXFEpCsHcXzWoqwS\nZjuGRTxKcyle6qDE+Xgb+8VRSnc/Z+BMvNlBxmSzUQTz65LxVNzkn5LsQBlhn2vNnOExEXAnYkQY\neX6y4erU3Jj4xU4ro0gJGopfBrNFHnUv7sQicYn9Luf2zTuP27QG9GnC+7KCYJpmwjaOX/Fiau4C\n4fKfdZDRiUdwesbzhSIY78+YW4ePUuOfcXvbO3PxjpL/VbpceMTsZDyGG0qsf1nkK2newDNllChJ\nQ/fpxXHCQNfnvZDXZnk9WbgiGa/S6h2dWCTykibDIjfrh+7AeuwWl1OTa/GtMHYmeYbaI87sJcl4\nFF+VUOZrkbM02STiW68DOdG52Is/k/Ey3CLCRW5omShn2I4rcbLyxep6PIT7E6XOEsbfUVLOZLBZ\n9MDuEjXrLBFq3u1W4Apx/O7BRRUUmyKC/NYKMorQUK0E6pph4Q27lauqH9eaaK4WN9Gp9amWyRaR\nnvSEV5XvZn4vShg4QfSc1tSpVC/Ia6c0mYcP/LsUmYhd4rpdKbzpVp0z+gEDBgwYUIG/Aayknkb9\n3TKcAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\sqrt{y^{2} + 9}$$" ], "text/plain": [ " ________\n", " ╱ 2 \n", "╲╱ y + 9 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z = sympy.sqrt(x ** 2 + y ** 2)\n", "z.subs(x, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "再替换 `y`:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAsAAAASCAYAAACNdSR1AAAABHNCSVQICAgIfAhkiAAAAMxJREFU\nKJHN0TFLQmEUBuDHFIKGFhVcmvwFDW35T/wHQXOLo3/AsSVozTlwCwcbGhoaxHLJTVAIAo0cbLjf\nhevHTW5bL5zh/d5zXt7vHP6AcsQn+MAyaOe4xhMW8fA2qm9cpGIlap5hgDrecIPxb7Ee9mU+2CfG\nKEX8EfeoYoMmrvCaNzzFSYa3MUejSKwyVugVjfoeamd1Q8kfWjnutdjhE885zl+Sy+7gTrKFLE4l\nl+zEzWfo4yjwEm4xwmH6kEULl1jjGC/ohij/BT8UwSaNbxctpgAAAABJRU5ErkJggg==\n", "text/latex": [ "$$5$$" ], "text/plain": [ "5" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z.subs(x, 3).subs(y, 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "还可以从 `sympy.abc` 中导入现成的符号变量:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAADgAAAAXCAYAAABefIz9AAAABHNCSVQICAgIfAhkiAAAAwFJREFU\nWIXt1luIV1UUBvCf5jRTpqNEmdeGSaickHoIqYx8CZoIwiAhUCh96KVMnHqIsIjCSzREQw9JEIE+\n2IUukJQPGjVRQSlRURFCklFZ2kOF0lDaw9qH/57tmWb+MzqCzAebc9Z31t57rX1Z3+Esx6RxnGsJ\nrkc7bsAT+GAc5z+tuACbM3sFjmLumQnn1GMxjuOyZE/HCZHoWYFJ4ohWV6JLJHjNGYuowDL8jvWn\naLxt6B3GpwUdQ3zrwFZswWuYhlmYOdqA7hIrvnW0A2RYg6f8f4Gbgsdwfs23DuzD7GT3oA+TsWmI\nPiPCwjTxWHCbSBDaDL1DPeLeljgXe7E64+7ED+l9Np7OO0xuIrj9+KcJ/xI3iWO0E5fgFo1dyDFN\nyMkXNd8eSH22Z1w75osj/bOo2AuaCawVl2Mprm6mY4ZO/CmOed6m1/jegQdr+DYcFvqZozeN1Zrs\nFXioLojl4m48jMfxXGpXYVca5KXC/8s06QbcjI14Fh/hutpUh0ef2O0Sd6cYugr+E/yR2Z14q+x8\nJfYU3EqDE/qssOFCsTM7Ne5WFeSBmiBHgtednAQR9DG8m7Xd+Bf9mV+7WHg07uBiXCzOf4U38Vdm\n5+8VjqTWiRcz/itciouGy6YGU/F3wZ0jdvUNcXer1ptyeC/zPSqKERoJ9osEfxQadT/Ow30jDOpz\ncXwqDGTBNovfnKxnc8XOfFzw3en5asbNxK+VUSX4k/jT2CFEvQ8Hhf6NBOWKjwXfYl7BzUrPrzNu\niigo/bIjKSrqN5VRJXitEN57k8MCvCyEvdX44h3cWHCVPP2Scd3i1G0ofJeKmoBGgl0G79ZBIabH\nMWNs8TaNvULrWjKuEvJch3vwAt4v+i8Ri4TBQr8OczJ7Hr7DoWS3qP+TqeNbimezeAarMvuIkJ4r\nkr1anKy1Rb9bRSEaKHgr8UgaeCOeFBo4X1SvT0URGRByMhW3i9Wu+A/Fyr8ifsxP4Hs8Osoku7Eo\nsxfhbTwvakRb4T8H94xyrglMYAKnCf8B3eKPw+iETXAAAAAASUVORK5CYII=\n", "text/latex": [ "$$\\sin^{2}{\\left (\\theta \\right )}$$" ], "text/plain": [ " 2 \n", "sin (θ)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sympy.abc import theta\n", "y = sympy.sin(theta) ** 2\n", "y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对 y 进行积分:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ0AAAAZCAYAAAArBywYAAAABHNCSVQICAgIfAhkiAAABZBJREFU\naIHt2nusXFUVx/FPoaWVlhQCCC0CV3ILVAIWGorWIqUJJkAaAgQrj5iC4h+mhEdjMEExQVM1iqaF\n8FCI4gMaFAIpFYRAqqVAoCUkGh5/+QAJCqEEDEjB1j/Wmdx9991z55w7U4fqfJPJzFmzH7+91tmP\ns2YYMGAX5qx+CxjQF/oW9yNxT786H9A3JhT33drYh3AzvotfY68O7ZyHOzLbAtyN+/EH3IqDmgrs\ngsV4HVf0oK0pwic5Q8p+OgD79KDfbmmnm7L2prp7FvchPI1Z1fVKrOlQ52l8KLk+Dg9i7+p6Bn6P\nf2jvhF5zLnYIx3bDZHwDe2b2Ie39tBu+Xajz36Sdbtprb6q7J3HfA1twUWI7B38dp+MF+HlmW4/h\nzHasuAnWjtNWrxkWzu+GlTgms9Xx0yx8v8u+u6Gkm87a6+ruWdy/gpcrYS2+WFWa0qbOapya2f4p\nBvHhzL4Vr7Xr/APIXrirYK/rp5twyE5T1552uqmnvY7uCcc9PdNNw5ViD96W2OcWyrbYHSfjocz+\nJ3E+mJ7Z3zV6Od5ZTMURWIR5XbRzCh7PbE389AiWddH/RCnppr72Trq7inu69XwO+xq7DH4Kb1UV\nc5ZgI97P7J8Qs+2VxDa7ErShOIzmnIlPilk0TWiHFZiDa/EZ3IblSZ1rxBayGk8I503H8WJLSoO1\n2NgVo4mfNuN8fC+xzcPVeAP/EoG+pSrb4gRcgr+LGO2NVXih5vhLuptoL+lO6Vnc78U7eCB5PYx/\nVx2U+EkluA7fqdpaWLP8eMwVszHlAvw0s20u2PYVDl6PLyT2NfhzVvZuHJXZmvhppniCa7EQb4qb\nosXP8Jfkeqm4ufZPbHMrbccm1+ONv6S7ifZcd05P4r67mHm3Z/bTxF5/TaGxaXgGk2p0PCz2+2/V\nFNqJZfij0amcGbg+K7fB2JuOCOBzRmv/khhrGuzfGn0wbuqnKUZWp0nV5/VZmR/izmQMr+GygubV\nwt90Hn+uu6n2VHdO13Fv7eEHibs7Pwe0Doq/KjR4usjF7OjQ8VQx0B/hazWE1mGjOKy+JJ6gLhFn\nhhUN2njGaO2tM056HnnV6LxVUz/tI9IFxJZ5uHhyTLkcn60+ny5W4ucLel/AxzFf5/HnuptqT3Xn\ndB331k13QPX+bPLdZOGMjcpL7bn4ZYeOJ4ml+AHlJO3RwglP1HzdVNV7WQRxrdiq1uDFSlNdSmfU\nnOfxkeS6qZ8OFisqI3mql8bp77DqPT8rwXvV+7DO4891N9We6s7pRdwxkktJzwFLK9tJhfIz8WSH\njoll9euZ7fM16nXieCMBIpx0mzgvTU3sG7TfXnP7cjHeocQ2Hz9Irpv66VKcUX3+dFXuykK5XMN5\nhe++Wn23WOfx57qbak91p/Qk7q2VrpUYTGfYSvwYvys0erY4rI7HhdiOb2b2RR3q1eEoo1e1F0Wy\nc7uRbHgv2CKedFv5q6Z+OkFsRfCYWOWWFMqdVfWzDm+LdE/O/Kr+ozqPP9fdVHuqO6Xncd9k5O6+\nSGx709o0/CAOHafjJeJA/IvstdbY3+omwnJxbpmd2A41dhZuqvrN+VvBfrGY9XMy+wKjs/d1/XSa\nyPSnnCLSJOkqsj9uTK7PFzfXrMT2UeHPk6vr5TqPP9ddV3tJd4uex/1juE+cm9YUxLQ4UHlWp2wV\nASy98hkwES7AVeLJb5VYzq8X2wyxXTxV9bdNpBemC4dvSeyPiuDeKf4csEMkOK/O+jtV+Id6fpot\nZnyJBfiNOGRfW+mfmZU5sfr+BlwnkrlHNxh/SXcd7ePp7mvcL8OXu2lgwC5JX+O+Cfv1q/MBfaNn\ncW/3f7p2DIsE4670o/2A7ulp3JvedKU/7Q3436evcV8nfm4Z8P/FIO4DBgwYMGDAB5v/ACgJsG1Y\n9yxvAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\frac{\\theta}{2} - \\frac{1}{2} \\sin{\\left (\\theta \\right )} \\cos{\\left (\\theta \\right )}$$" ], "text/plain": [ "θ sin(θ)⋅cos(θ)\n", "─ - ─────────────\n", "2 2 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Y = integrate(y)\n", "Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "计算 $Y(\\pi) - Y(0)$:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJMAAAASCAYAAABfCexoAAAABHNCSVQICAgIfAhkiAAABchJREFU\naIHt2WuMXVUVB/DfQEsZWoHwHAkNLS3iI1FKKkVAYkRISEggTeQRG0C+CBgCRh5KBIsFITXRiGCx\nEBiF8LAEmkAT1AhaPwA18hAioZKYajHIEFCgUrTt+GHt4+x7Zp/be+7cfiHzT07u3Wuvs9ZZ+6y9\nHvswjWkMCEMN9L3wFD7ZQtbLWIEnsAWL8Q1cgpcSz6m4Dc/jPbyPHZmM3+PW9P9I3IhNGMcBuAL/\nyPjbyDsc38ZWbMdsXIXXCrZ8BMszme+l8TsZz7H4SpobFmv23fQsOY5J6zCMQ7EB1+LVPvTuKntz\nXI6ZYu1z9CXv0+mhxneitI7x2vUffLXG8/UCX36dmvj2wWYsy+69Gi9ijz7kzccbWJrdu0y8lBm1\nZ/w4/orj0ngEfxEvtcIiPIo9M9oqvI2jMtrR+CX2TeM5WI/XMa8PvbvC3hyHiUCwvEZvLe9jWIdR\nEZXaOtMmrMbDWJnk1bEKc4Xn75bRj8ct2fgGseD5g+6H/+KiPuStxZs1nmGxy87PaDNEhL0so83F\nGC7NaD8Q63NWRjst0W7OaOuwUCcWJb77+9A7aHvrWJ2ebXmN3q88hEO1dabf9MBzS4E2B4+JNFFh\nIx4p8L6Ax1vK20M44YYC78siwlS4QETUfQu8Oc7Fv3ByRjtbrNlNGe1dEW0Oqt3/ltjpbfUO2t4c\nS03YsDyj9yvv/xi1a5yphFVYko0/lHTfWuD9hXiJbeSNJHnrC7wbRASs8Cv8aSfym7AS2/CpjPaC\nqGvm13hfE+lkEHqnYm+FObgj/a87U8/yuuXPtpiFa7C/8OQFovDc2OWe47E7ns5oh6Xftwv8W7B3\n0vV+j/LG8G+d9U2FQ3CgWIftOEEs0Ik4RSzyPFyHZ7vYMR/n4WKdBfixYnPkReohONjE5huagt6p\n2Lsto3/T5IJ7KvI6MKp9ZHpF5PMKy0TnNdLlnj+YvGuPS7qvK/DfneYObiEP7hRpJe9eR0QnNC7S\n0AHp/0u4MOP7nOimPlGQe5roYP8oGoTdCjx13CQctyq0+9FbYSr2VjhKOFOFUs3URt4kjGrvTPWF\n3F149M0FXjhJObQvUTaIKFrHxY7oVR6xezaKVp7YSdfjmSRvf+Gg46KoHK7dv1n32mCGSFVPCedo\nwkJRR12f0frVO1V7iXf2U50dcmnte5LXy07qFTtq4+0iRJ7ewH+xWPw6xrromJ1+3ynMNcmrZC4R\nu+n7YrFuF8+8VXQqbybeV8QZTo7Notie1SB/m4hQS8Q5UAmzcK/omL6V0fvVO1V7CecYFcV/N/Qq\nr4hR7SLTevyuQN+clNUxU+zQGwpzs9ND/rAw97johNrI64a/4dfZ+HVlO34r1uPDafxRkw909048\nO0TNk2NIONJ3Gp6jV70VBmHvCH5U4GnKCjuTN7DItMjkBSRC/qYC/RjhNG8U5raIonNuYW4hnmsp\nrwkHihPpNRmtKU1VBf+YcJpnRYhfkPFsT79DIsXnWCFqomsz2rkt9eYYhL1fEJtibXatS3Nnp/FS\nzSitXxGjukemI3Tm9zUm8nCF6nDumsL956W5iwpzxA5+VWfBtyDdUz9V70XeZWL3H5rRLheRM+9S\nzhF1Xm7bEP6Jn6fxnqJb/bM4SK2wOD1D3lnBl5Uj0uqWenMMyt465ilHpp7kNUWmimGvwtyJ4rBq\nbUZbiZ9k/EP4Gp5Mc3VU1X9Trl4ldl7+OeUSUXDe3oe8OeJlVe3rIvGd70ydafgB0drnDvtF4TxX\npvFWfE80FnnKvVSknvwFfz7xHo57sut+cWTQRm+OQdlbx8zab9/yDhKHgi+a+NYzJj7cfinjO1K0\n/D+u3f9ZEaF+JhxtheZdcIbYdYsbzYqWdZ34dHEHHlJOfb3IGxYv9S48KDqvzzTw7ici8xpR59yn\nM51VOD/N3Z3kPSA+1OZ4S/N3tBV96mWw9hKp+wn83URn+WTS04+8aUxjGtOYxjQ+aPgf1hMJbOhz\nlqIAAAAASUVORK5CYII=\n", "text/latex": [ "$$1.5707963267949$$" ], "text/plain": [ "1.57079632679490" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "np.set_printoptions(precision=3)\n", "\n", "Y.subs(theta, np.pi) - Y.subs(theta, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "计算 $\\int_0^\\pi y d\\theta$ :" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAACEAAAAZCAYAAAC/zUevAAAABHNCSVQICAgIfAhkiAAAAc5JREFU\nSInt1k2IjlEUB/DfMOWzMQtponzUu5gtZVZYCCVlYTaMmrKwsdDEUlamLJSyYEXKR5RCDRELUT6i\nkFnZDQtJioWJmfKxOM+bd67rfXg9Q03+9dS5/+fcc//3nOee5zJFseVfC+jG5d+dNK1iEX04l3A9\nuIhrGMYJLKp43Ql4jFkN4xW4gc5iPBd38AZLJ0NAD04n3FXUEm45vuL8ZIg4go0J9wEvsSDh3+Ft\n1QKm4xnaE34YY1iW8K8xmgvUIVL0UaQr93zB6szc9Tia4eegK+EWFrFu1Ym68jacwhOcFGf9EUYw\ngGMYxyfcyyzWh+MZftSPO95dbGZf6txvYj0viBTDUCZ4I2biabGRMtTEdzJY5tiJS4U9Gw9K/Htx\n8BcEzMBDHE5f5JrVdtwv7G58Lgm+DWdLfNpEma9jT4kvIrUrC3sTnjfxnSd2V4ZB7E+4/rqRZmIN\nlojOR5yYxaLuOfSKltwMO8SHeCDhV9WN9FzvxU3fS/CqELAOVzILbMXOJgLW4pAow5kGvl0c0yxG\nsKFh3IEX2Jzx7cLtJgKIzviznpNmpiUMYFcVgf4EdzG/ikCt3idqeK+in1CrInKXl7+OIXFB+Y+p\nh284OVxxrb4HkAAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\frac{\\pi}{2}$$" ], "text/plain": [ "π\n", "─\n", "2" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integrate(y, (theta, 0, sympy.pi))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "显示的是字符表达式,查看具体数值可以使用 `evalf()` 方法,或者传入 `numpy.pi`,而不是 `sympy.pi` :" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJMAAAASCAYAAABfCexoAAAABHNCSVQICAgIfAhkiAAABchJREFU\naIHt2WuMXVUVB/DfQEsZWoHwHAkNLS3iI1FKKkVAYkRISEggTeQRG0C+CBgCRh5KBIsFITXRiGCx\nEBiF8LAEmkAT1AhaPwA18hAioZKYajHIEFCgUrTt+GHt4+x7Zp/be+7cfiHzT07u3Wuvs9ZZ+6y9\nHvswjWkMCEMN9L3wFD7ZQtbLWIEnsAWL8Q1cgpcSz6m4Dc/jPbyPHZmM3+PW9P9I3IhNGMcBuAL/\nyPjbyDsc38ZWbMdsXIXXCrZ8BMszme+l8TsZz7H4SpobFmv23fQsOY5J6zCMQ7EB1+LVPvTuKntz\nXI6ZYu1z9CXv0+mhxneitI7x2vUffLXG8/UCX36dmvj2wWYsy+69Gi9ijz7kzccbWJrdu0y8lBm1\nZ/w4/orj0ngEfxEvtcIiPIo9M9oqvI2jMtrR+CX2TeM5WI/XMa8PvbvC3hyHiUCwvEZvLe9jWIdR\nEZXaOtMmrMbDWJnk1bEKc4Xn75bRj8ct2fgGseD5g+6H/+KiPuStxZs1nmGxy87PaDNEhL0so83F\nGC7NaD8Q63NWRjst0W7OaOuwUCcWJb77+9A7aHvrWJ2ebXmN3q88hEO1dabf9MBzS4E2B4+JNFFh\nIx4p8L6Ax1vK20M44YYC78siwlS4QETUfQu8Oc7Fv3ByRjtbrNlNGe1dEW0Oqt3/ltjpbfUO2t4c\nS03YsDyj9yvv/xi1a5yphFVYko0/lHTfWuD9hXiJbeSNJHnrC7wbRASs8Cv8aSfym7AS2/CpjPaC\nqGvm13hfE+lkEHqnYm+FObgj/a87U8/yuuXPtpiFa7C/8OQFovDc2OWe47E7ns5oh6Xftwv8W7B3\n0vV+j/LG8G+d9U2FQ3CgWIftOEEs0Ik4RSzyPFyHZ7vYMR/n4WKdBfixYnPkReohONjE5huagt6p\n2Lsto3/T5IJ7KvI6MKp9ZHpF5PMKy0TnNdLlnj+YvGuPS7qvK/DfneYObiEP7hRpJe9eR0QnNC7S\n0AHp/0u4MOP7nOimPlGQe5roYP8oGoTdCjx13CQctyq0+9FbYSr2VjhKOFOFUs3URt4kjGrvTPWF\n3F149M0FXjhJObQvUTaIKFrHxY7oVR6xezaKVp7YSdfjmSRvf+Gg46KoHK7dv1n32mCGSFVPCedo\nwkJRR12f0frVO1V7iXf2U50dcmnte5LXy07qFTtq4+0iRJ7ewH+xWPw6xrromJ1+3ynMNcmrZC4R\nu+n7YrFuF8+8VXQqbybeV8QZTo7Notie1SB/m4hQS8Q5UAmzcK/omL6V0fvVO1V7CecYFcV/N/Qq\nr4hR7SLTevyuQN+clNUxU+zQGwpzs9ND/rAw97johNrI64a/4dfZ+HVlO34r1uPDafxRkw909048\nO0TNk2NIONJ3Gp6jV70VBmHvCH5U4GnKCjuTN7DItMjkBSRC/qYC/RjhNG8U5raIonNuYW4hnmsp\nrwkHihPpNRmtKU1VBf+YcJpnRYhfkPFsT79DIsXnWCFqomsz2rkt9eYYhL1fEJtibXatS3Nnp/FS\nzSitXxGjukemI3Tm9zUm8nCF6nDumsL956W5iwpzxA5+VWfBtyDdUz9V70XeZWL3H5rRLheRM+9S\nzhF1Xm7bEP6Jn6fxnqJb/bM4SK2wOD1D3lnBl5Uj0uqWenMMyt465ilHpp7kNUWmimGvwtyJ4rBq\nbUZbiZ9k/EP4Gp5Mc3VU1X9Trl4ldl7+OeUSUXDe3oe8OeJlVe3rIvGd70ydafgB0drnDvtF4TxX\npvFWfE80FnnKvVSknvwFfz7xHo57sut+cWTQRm+OQdlbx8zab9/yDhKHgi+a+NYzJj7cfinjO1K0\n/D+u3f9ZEaF+JhxtheZdcIbYdYsbzYqWdZ34dHEHHlJOfb3IGxYv9S48KDqvzzTw7ici8xpR59yn\nM51VOD/N3Z3kPSA+1OZ4S/N3tBV96mWw9hKp+wn83URn+WTS04+8aUxjGtOYxjQ+aPgf1hMJbOhz\nlqIAAAAASUVORK5CYII=\n", "text/latex": [ "$$1.5707963267949$$" ], "text/plain": [ "1.57079632679490" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integrate(y, (theta, 0, sympy.pi)).evalf()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJMAAAASCAYAAABfCexoAAAABHNCSVQICAgIfAhkiAAABchJREFU\naIHt2WuMXVUVB/DfQEsZWoHwHAkNLS3iI1FKKkVAYkRISEggTeQRG0C+CBgCRh5KBIsFITXRiGCx\nEBiF8LAEmkAT1AhaPwA18hAioZKYajHIEFCgUrTt+GHt4+x7Zp/be+7cfiHzT07u3Wuvs9ZZ+6y9\nHvswjWkMCEMN9L3wFD7ZQtbLWIEnsAWL8Q1cgpcSz6m4Dc/jPbyPHZmM3+PW9P9I3IhNGMcBuAL/\nyPjbyDsc38ZWbMdsXIXXCrZ8BMszme+l8TsZz7H4SpobFmv23fQsOY5J6zCMQ7EB1+LVPvTuKntz\nXI6ZYu1z9CXv0+mhxneitI7x2vUffLXG8/UCX36dmvj2wWYsy+69Gi9ijz7kzccbWJrdu0y8lBm1\nZ/w4/orj0ngEfxEvtcIiPIo9M9oqvI2jMtrR+CX2TeM5WI/XMa8PvbvC3hyHiUCwvEZvLe9jWIdR\nEZXaOtMmrMbDWJnk1bEKc4Xn75bRj8ct2fgGseD5g+6H/+KiPuStxZs1nmGxy87PaDNEhL0so83F\nGC7NaD8Q63NWRjst0W7OaOuwUCcWJb77+9A7aHvrWJ2ebXmN3q88hEO1dabf9MBzS4E2B4+JNFFh\nIx4p8L6Ax1vK20M44YYC78siwlS4QETUfQu8Oc7Fv3ByRjtbrNlNGe1dEW0Oqt3/ltjpbfUO2t4c\nS03YsDyj9yvv/xi1a5yphFVYko0/lHTfWuD9hXiJbeSNJHnrC7wbRASs8Cv8aSfym7AS2/CpjPaC\nqGvm13hfE+lkEHqnYm+FObgj/a87U8/yuuXPtpiFa7C/8OQFovDc2OWe47E7ns5oh6Xftwv8W7B3\n0vV+j/LG8G+d9U2FQ3CgWIftOEEs0Ik4RSzyPFyHZ7vYMR/n4WKdBfixYnPkReohONjE5huagt6p\n2Lsto3/T5IJ7KvI6MKp9ZHpF5PMKy0TnNdLlnj+YvGuPS7qvK/DfneYObiEP7hRpJe9eR0QnNC7S\n0AHp/0u4MOP7nOimPlGQe5roYP8oGoTdCjx13CQctyq0+9FbYSr2VjhKOFOFUs3URt4kjGrvTPWF\n3F149M0FXjhJObQvUTaIKFrHxY7oVR6xezaKVp7YSdfjmSRvf+Gg46KoHK7dv1n32mCGSFVPCedo\nwkJRR12f0frVO1V7iXf2U50dcmnte5LXy07qFTtq4+0iRJ7ewH+xWPw6xrromJ1+3ynMNcmrZC4R\nu+n7YrFuF8+8VXQqbybeV8QZTo7Notie1SB/m4hQS8Q5UAmzcK/omL6V0fvVO1V7CecYFcV/N/Qq\nr4hR7SLTevyuQN+clNUxU+zQGwpzs9ND/rAw97johNrI64a/4dfZ+HVlO34r1uPDafxRkw909048\nO0TNk2NIONJ3Gp6jV70VBmHvCH5U4GnKCjuTN7DItMjkBSRC/qYC/RjhNG8U5raIonNuYW4hnmsp\nrwkHihPpNRmtKU1VBf+YcJpnRYhfkPFsT79DIsXnWCFqomsz2rkt9eYYhL1fEJtibXatS3Nnp/FS\nzSitXxGjukemI3Tm9zUm8nCF6nDumsL956W5iwpzxA5+VWfBtyDdUz9V70XeZWL3H5rRLheRM+9S\nzhF1Xm7bEP6Jn6fxnqJb/bM4SK2wOD1D3lnBl5Uj0uqWenMMyt465ilHpp7kNUWmimGvwtyJ4rBq\nbUZbiZ9k/EP4Gp5Mc3VU1X9Trl4ldl7+OeUSUXDe3oe8OeJlVe3rIvGd70ydafgB0drnDvtF4TxX\npvFWfE80FnnKvVSknvwFfz7xHo57sut+cWTQRm+OQdlbx8zab9/yDhKHgi+a+NYzJj7cfinjO1K0\n/D+u3f9ZEaF+JhxtheZdcIbYdYsbzYqWdZ34dHEHHlJOfb3IGxYv9S48KDqvzzTw7ici8xpR59yn\nM51VOD/N3Z3kPSA+1OZ4S/N3tBV96mWw9hKp+wn83URn+WTS04+8aUxjGtOYxjQ+aPgf1hMJbOhz\nlqIAAAAASUVORK5CYII=\n", "text/latex": [ "$$1.5707963267949$$" ], "text/plain": [ "1.57079632679490" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integrate(y, (theta, 0, np.pi))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "根据牛顿莱布尼兹公式,这两个数值应该相等。\n", "\n", "产生不定积分对象:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAF8AAAAxCAYAAAC8snXfAAAABHNCSVQICAgIfAhkiAAABXZJREFU\neJztm2tsVEUUx3+1LaW20OJbaLGRh1qjSHw0sShtPxAhRqOJjcQX2kS/IBirn4ya4AsfQGjUSPzg\nE1MVNSS+FalSRbQYYsVX8REh8QVobERbEPzwn5udvd3unXt3t5vN3l/S3Duzc+ecPXtm5pyZW4gp\naJYAe4HPgJPyrEtR0QocAM4FBoDe/KpTXHwCvA5UAt8B3flVp3hoAg4B1+RbkWLkYWT8Y/OtSDHy\nE7A930oUModFfG4GUE+8wGZEVOO3meuWbClSjEQ1fqu5bs2WIjHu/Az8C5TnW5FiYyaKcvryrUih\nUxbhmfPNdVs2FfHRhLLmGqAZuAv4IIfyCoankecvyVH/1cByq9wO7AOm5EheQTGAjN8a1DAipwMH\ngWmmPNHIa8+RvILhaGSIQ+Y+F5SgaafElE818mbnSF7BcBEyxK9p2rSgLeabsyTzGWCFQ7tyoCFF\nfQOwBrgfWAdMQFsik7KjHitRtj/abJBKfiSWGyHvpmmz0LRZE1WIRQfwAIlRMBplwJ3A4b76BnTO\ncLwpdwJdKL+5L0X7qCwGhtDurov8SPQgw64MaDedaJGUzYXI+ADjSe3VHp1orbAZh5LA66y6y5CX\nggzyUIY6enQzcqslSH6oDLcUOMvcfx7Qdgc6ZInKXDQ1vAYcB1xAwnv8TEBhqV+npeaZZ626GrQn\nVY4SxWpgagZ6epwHvB9Sfii8hS/d4leBjhLnAGeEFWA4ERi0ZHl/E0dpfylwi69uPLAb5Qc2K0xf\nFabcDtwaUU+PaabPeWHlh5kazjbXAyjLXQj8YQQdaT57zAiYBzwFLDL1lwDLkCesBj5Gi1OV6bcT\n2Gzafk+4RakFeMlXd7nRyX+y1ox+2CFT7gOuAB4MIa8NuAH4EagF+pFNPowg35lH0K82BLzn++xK\n4Emr3OcrY5QZRFNJh1Xfhb5IVF5Go9JmPfAP8Kb1twH4D9hktatBxnOlA0V6daY8FSWAm33tXOU7\nswUZ/wfgC5K9sxqdbHn0MNL4ICN/RXL0cj2Z5Q1voQXeoxT4E3jO126BkbPMqisHvnGUMwvYj0a8\nzW6SM3Jn+a4Lbilwmrl/AjgG2IVi8BtRiLXYsa9tRgmPYXOtcnzez+8kx+xTkEf7vXG+ub5o1U0C\nfnOUcw8aueusukY0ou3F1lm+q/FPJhHDvo0y0G4033YBOxnpEaMRer4L4GsS0wAkzpS/tOrK0OK6\nieRpph6NxCBqUcT1DvJ+jxY0ldhhprN8V+N7kct+NFRL0KJTj+a951FSVZHy6dzyBgr1PLwQ9xer\nbj4arbf7np2D1qAgpqPR7/fmVpREDaIoLZR8V+PPMtd+FFrZXr4TJRIHkYeMNVtRFOXFzl4SY+cZ\nncDjjIzFm9CPF8Rfvr5BM8Fcq8+lYeW7Gt+b770z25uAydbndcC3JPZ8ykkdxqaqL/ddo7AKuMrc\n7wE+QlMlyDEqGLkFvgB4hcSak44BlMQ1WLo+avrdgbza++6u8p3ZhRbJRSisvA194XuBu1GkU488\n4VPTdhiFpFXAxchDvfpe5K0voE04L4q6I6qCaGg3mvtG4FWUd3ShXMRmMnBtyP5noilqlelzBrLH\nRuTVR1htg+Q7U0siyzwlaicx0WhGht9D8O5iTAhc5nwve+wlOT6PyZAwxo8PsLOMi/G9MDN+NTAP\n7AX+Jn5Basw5Ac3z/l3MmCwQNO2caa4bc61IMRJk/HPMdUOuFYkZSQ/a18j0MDwmJJXoNMZ/RBeT\nJdJNO21oT2L9GOlS1KxGO3fjTHktOiDP1otFMWnYjv6XthLtUA6R2S5jTACl1v1s9EpHP3rtYxi4\nGh2TxeSYo9D57D70X+V16ZvHxMTExBQY/wOjDi/8nWy6FQAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\int \\sin^{2}{\\left (\\theta \\right )}\\, d\\theta$$" ], "text/plain": [ "⌠ \n", "⎮ 2 \n", "⎮ sin (θ) dθ\n", "⌡ " ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Y_indef = sympy.Integral(y)\n", "Y_indef" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "print type(Y_indef)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "定积分:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHIAAAA4CAYAAAAhHqx3AAAABHNCSVQICAgIfAhkiAAABoBJREFU\neJztnGmMFEUUx3+wF8gtXhyLK5cCKhpZiYLAkmiAGIgmokSIKInEBIG4+oEYNcELDyDgEQkfPBCD\nihoSiXggKCgSQVFAUUBBiIACHnhwyfrhX5Xp6e2d6e6Z2R5C/ZJNzdRU13vT1VX16r03Cw5HAzQB\nngI2AfuALcBeUy5MUC9HFiYDB4EvgPOBsUAfoByYCjQFpiWmnSMUNcBx4EpgK7Da89lIYBAa1EmN\nr9qpRdMcr38MeA/40vS12/PZOOAzoCtabh0FpDSHa/sD1cB44F+gm+ez3kAz4ChwGnBRDnIcIchl\nRo4z5bKAz8YAr5jXm9De6ShSfgI2J62EQ8SdkT2AStKNG0eCxB3IoaZcmy9FHLkRdyBrTLk+X4o4\nkmEPcBgoS1oRR3x6AnXAuqQVcaSIc44cZMoN+VSkwPRH3qc2wADgQeDjRDUqAl5CM3Jy0oqEpCUw\nw/N+NPAP0CkZdYqHrWgga7I1LBIuBk6Q8jy1RvqPTkyjIuBMdBPqzOuTgSZoabX+3j5I/0sT06gI\nGIluwr5GlDkEhcnuylN/C4CZWdqUAVUB9VXAPBQsWAy0MvVnA+3yox6zkNcsaNVrSH5kZhgBH8Tt\nIAZjjMx5eehrAvA4maMxpcADyNnvpQrFXDuY97XAXPO6KfBowDVxmQQcAZqHlB+ZleimzorbQUy6\nk1ukBuBaNJCgyExVA+1q0b7qpRw5P27z1N2AZo6lA/BkjjpaFpHu/swqP4pnpwToZ15/HVPBuGxD\nAey4DEbL31LgHGAYqSfbSyt0TPF/vymm/cueujbI32ydInuQhdwlBz0tVwEfRZQfGmskNKahUIFC\nYAOBS2L20RU4REp3+9c6oO31wN2+umbAfnT29DLT9FPhqRsN3BNTT0s30+81UeRHWa6qTXkc+Ca2\nmvW5DrgC+A0p3d7UT0JRlpnoS72Igtj2munoKZ2DMhFqgBZGz1pgjWn7A+GNgiHAG766m4xOi3z1\nA9ADcsRTtw64GXgipDxQAGIisANoC2xE9/iTGPJD8Qx6Ar6KemEGegEf+urGAi/46tYF1LVHX2Qp\nqb0PZADsiKnPm2jl8bIEZUAs8/wtB/4DVvnatkEDEZYJ6ATQ2bzvgpwVazxtosgPxVo0kM/HubgB\nbkQZBN4Z0xJ42tduJfUHEjRg35Juhd5O/HPuu8iwspQAv5PKdrCMMDKm++rLgO9CyuoLHENWuZf9\npDxRoeWHNXZKSOXd5DN0tQo4CyVtLQDuRCZ3lKy7DehLWY6askUMfX4l/TzYCc2yNb52w035uq++\nHfBLSFkPoxVlsaeuN1pprKETWn7YgbyA1Jkmn1GPn5HXZRHan+YCu6j/lGYi8h6RgS2kljmQpQvp\nNkEpMmpWUX8ZrUQrRDbaIsv5fTQrLUPQkmmPHqHlhx1IazEeI79Rj2q0LE5EN6EL8Co6/FdkuK5Q\nvINMf4s98uz11A1Hq8h9AdcPRHt2NrqjVc4/02rQof8QsrZDyw87kH1NuREFlPNFH9Jn3y506D2B\nntrGZj2yhO3ZzB64vWfYWmA+6ec8S3/0MGTjT1//oBVvsKffKVHkhx1Iuz8WIkdnKtDR874z8D3p\n/twygo9KQfVlvjIqs0mleh4APkVbC+ghqyA4hDcCeIvUHp2JrcjpUOXR9VnT9zY04/ZFlB+K3cig\nGB+3gwYYC9yLbt4jwEPIYq00nw8GPjeyj6KjSgtgFJo9tn41mkmvIQd7HfAjcH9MvYYjwwNTvg08\nh/bwZgHtOwK3RpTREy3Ds02/PdD9XYFm3OkR5IeiLSlvSK+4nTiSZwAaxAO433AULWH2SOvpWE36\nec1RREQZyFM+WamYCTOQ9ujhfh5wknMQ+BuXjHxScy7aF/0RCkeRkW1pvcyUKwqtiCM3sg3k5aZc\nXmhFHIVlJfIL5pr45EiQ5igy7U99cBQhmZbWocivt6SRdHHkiTnII19u3i9ECVFBSbf9kAP3FuTg\n7R7QxpEQm4HtaEmtRJH3oOhBObCTVF5oNYpQOIqE+SiM1AH9JGATwVH6q0n/bx5NgL+A8wqtoKNh\nvHvkNHTc2I5ifMMIzoepQpEQSx1agi8sjIqOMHiPFftJZTdn4gyUe+nlMDn8MsiRO3H+q8cf1I9L\ntkQPgiMh4gzkFvRDGEspyufcmReNHI1GKcpHtb86sil8jgQpiXHNCZTbegfKih6FMuEOZLrI4XA4\nHA6HIwn+B7jjVvTlLJflAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\int_{0}^{\\pi} \\sin^{2}{\\left (\\theta \\right )}\\, d\\theta$$" ], "text/plain": [ "π \n", "⌠ \n", "⎮ 2 \n", "⎮ sin (θ) dθ\n", "⌡ \n", "0 " ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Y_def = sympy.Integral(y, (theta, 0, sympy.pi))\n", "Y_def" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "产生函数 $Y(x) = \\int_0^x sin^2(\\theta) d\\theta$,并将其向量化:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Y_raw = lambda x: integrate(y, (theta, 0, x))\n", "Y = np.vectorize(Y_raw)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEVCAYAAAD91W7rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4U+W1x/HvYlBUFFQUFVDqgMOtKFARxdagSAEtVkUF\nZ7FXbatWqx2w3Epbr1VbxaF1QrAgjq0IDogWL1FBZZBBEFGEqjgAUsSCiAJn3T/egDEcTpJzkrOz\nk9/nefKcDG92Fvs5rLxn7Xcwd0dEROKrQdQBiIhI3SiRi4jEnBK5iEjMKZGLiMScErmISMwpkYuI\nxJwSuYhIzCmRi4jEXKOoAxApBjNrBiSBVcBxwM7AgcCxwFJgnrv/M7IARQpIPXIpV78CngeaA/sC\nuPvzhIR+F/B/0YUmUlhK5FJ2zKwBcB7wN+B77v4GsNrMWgLLgK2B7SILUKTAlMilHB0ObHD3ue6+\nMvXcICABfA4c4+7/iSo4kUJTjVzKUTdgcvoT7v7LiGIRKTr1yKUcHQ1MizoIkfqiRC5lxcwaAkcA\nM6OORaS+KJFLQZnZt3Jos7uZbVukEA4FmgKzi3T8GpnZfmZ2kpldY2Yds7Td4rkq8jmSMqNELgVj\nZnsDXXJo+glQrJr1EcBSd/93kY6fzQnAh8DNwFVbapTDuSrmOZIyo0QuhXSxuz+UrZG7rweeNrNz\nihDD4cDcIhw3J+4+xN2nAm2Af9XQ9KKazlWRz5GUGSVyyZmZfcfMnjCzj82sZ9rzB5vZp8AFuZYD\n3H0a0L0IYXahCInczK43sx55vOUk4H+3cKxDgA8ynutqZgPM7Coz2xFqPkdm1s3MPjKzNjUdQyqD\nErnkzN2nA5cCOwJT0176L+DPhMk3a/I45Cdmtm+h4jOznYB9gDmFOuZG7v5rd38uxzj6ALcBrbbQ\n5ATSZpamzsF57j4ceA/om9Z2S+foJUIJaXEOx5Ayp3Hkkhd3f8/MXgHOBG43s17A20A/4I95Hm42\n0Al4Z+MTqdrxf9fwnlfdfewWXuuU+lnwRJ4rMzsJuJrwhZek+l75YcB1aY9vSHu8P7Ah7bXNzlFK\nJ745xLKmY0iZUyKX2hgBXGJm04Gv3H2GmW3r7lXpjVI90w3AdwnJtSfwv+4+P9XkU6Bd+nvcfREw\nsJZxdQKqqENpxcx2AM4BFgEHEKb59wBOdvfTzKwToQ6/BzAdaAgc7+4DUvE/Djye5WO2dXdPfd4e\nhMTeycy+Q/iC/ENa203nyMzaAWcTJjtdDDyU4zGkzKm0IrXxGHAQ8K3UQlQQEtomZrYnYYXBpwmr\nDz4NPAK8n9bsC2CrAsbVAXjb3b+owzF+SLhQ+TLh39gBGA9sHCq4CzAf+C93HwOMJkxAykf6uToG\neMrd7wFGEcox6SWcL4CtzGw74FHgJncfD+zA1+WZbMeQMqceudRGY2CNuz+Y9tz69Abu/j5AaqGq\nVak1T57KOE4zYEX6E3UsrRxC3Wd0PkNIzHMIfz08b2aXE3rmuPt4M/sjcH+q/RG1+Mz0c9UaeDN1\nvw/wtLsvT3t94zk6GZjj7ivNrAnQ1N0/yfEYUuaUyKU2jgYmZTy3xMyauvtqADM7gLDKYEfgxdRz\nJ7h7ejLfna8TEFD70oqZbUNYrvaufN+bdozDgQvc/QIz25VwwfIuoD/Qw8yOT/2F0Q24PvW2c4Ch\nZtYz1VPORfq5+iR8tBmhbHJxRtvdCX8BHMzXk5yOA15N+8xsx5Ayp9KK5MXMTgcuAxqbWee0l14A\n0h/3IIzOMKBJ6iLgsozDHUrG4lZ1cCDh97kuPfJlwGup2v4ZwJWp5xcR/i1TUsMrV7r7Z6nXPgd2\nJeMviyzSz9XfgfbAj4Dfb/xLJs2hhC/Nh4DWqYvLLYB1wPY5HkPKnKWuuYjUiZk1B65y90E5tm8C\nXOfuPy/Q558NDAd2qGONvOhyPVeFPkdSvmrskZtZEzObYmazzGxeqjaY2SZhZp+Z2czULaf/yFJe\nUjXw5WbWIse39APuLmAIBwKzSz2JQ17nqtDnSMpUjYnc3dcC3dz9UMKfbt3M7Khqmr7g7h1St2uL\nEajEwq2EGY01Ss1G/NTd3yrgZx9MmCQTFzWeqyKdIylTWS92ps3U24owbKq6WqAVMiiJp9TY6KE5\ntFsMLC7wxx9MGN8eC9nOVZHOkZSprBc7zayBmc0i7Dw+0d3nZTRx4Egzm21m48zsoGIEKrIlqan5\ne1K4C6cisZI1kbt7Vaq00hr4npklMprMANq4+yHA7cCYgkcpUrPvAG+5+8dRByIShbxGrZjZ/wBf\nuPufa2jzL6CTu2dO9NDwGBGRWnD3GsvX2UattEgNldo44eI4MrbQMrOWqYkIpMYVW2YSTwsmtrdr\nrrkm8hgqMfZc4j/ggAOYOHFi5HFW6vkv9Vvc489FtouduwMjzKxBKunf72HK8kWpxHw3YbnMH5vZ\nemANYciUSFFNnDiRe+65h+7du9O0aVMSicSm11asWMHQoUPZddddad++PZ06ddrygUTKQI2J3N3n\nEKZYZz5/d9r9vwJ/LXxoIlvWrl07pk+fzpw5cxg9evQ3XhsxYgTdunWjY8eOnHvuuTzwwAMRRSlS\nP7TWSo7Se3xxE+fYofr4W7VqxYIFC6ptv2jRIvr27UujRo1YsSKfmfPFUY7nP07iHn8utNZKjuL8\nyxDn2CH/+KuqqmjYMKwUm7p8E6lKO/+lJu7x50KJXMrO/vvvz9KlS1m7di077LBD1OGIFF29LZpl\nZl5fnyWV7d///jfDhw+nWbNmHHzwwRxxxBFRhyRSa2aGZxl+qEQuIlLCcknkKq2IiMScErmISMwp\nkYuIxJwSuYhIzCmRi4jEnBK5iEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjMaWMJ\nEalY7lBV9fVt42OAbbeNNrZ8KJGLSNlwh+XLYeFCWLQo3BYuhCVLYNWqzW9ffgkNGoSb2df3DzgA\nZsyI+l+TOy1jKyKx9fnnMHkyTJwIySTMnQtbbw177/3N2x57wPbbb35r0iQk8FKm9chFpKy4h8Q9\nfnxI3rNnQ8eOkEhAt27QoQM0bx51lIWlRC4iZWH5chg5Eu65J5Q+Tj45JO8jj4xXLbs2cknkqpGL\nSElyD+WSoUNh3Dg48US4917o2rX0yyH1rcYeuZk1AV4Atga2Asa6+8Bq2t0G9ALWAOe5+8xq2qhH\nLiJZucOTT8Kvfw0NG8KFF8JZZ8GOO0YdWTTq3CN397Vm1s3d15hZI2CSmR3l7pPSPqQ3sK+772dm\nhwN3Al0K8Q8QkcoycyZceSUsXQo33QQ9e6r3nYusE4LcfU3q7lZAQ2BFRpM+wIhU2ylAczNrWcgg\nRaS8ffghnH8+9OoFp50WLmL26qUknqusidzMGpjZLGApMNHd52U0aQUsTnv8AdC6cCGKSLn66iv4\n3e+gfXvYbTd4+224+GJopKt3ecl6uty9CjjUzJoBz5pZwt2TGc0yvzerLYYPHjx40/1EIkEikcgn\nVhEpIwsXQv/+sMsu8Npr0LZt1BGVhmQySTKZzOs9eQ0/NLP/Ab5w9z+nPXcXkHT3h1OP5wNHu/vS\njPfqYqeIAPDgg/Czn8GgQXDZZSqh1KTOFzvNrAWw3t1Xmtk2wHHA7zKaPQFcAjxsZl2AlZlJXEQE\nYPVquPRSePllePbZMJlH6i5baWV3YISZNSDU0+939+fN7CIAd7/b3ceZWW8zewf4HDi/uCGLSBzN\nmgX9+kGXLqGU0rRp1BGVD83sFJGie/RR+OlP4ZZb4Mwzo44mXjSzU0Qid9ttcOONMGECHHJI1NGU\nJyVyESkKdxg4EMaMgUmTNCqlmJTIRaTg1q2DCy6ABQvCaoU77xx1ROVNiVxECmrVKujbN6wL/vzz\n5b86YSnQnp0iUjDLl4d1wffcE0aPVhKvL0rkIlIQK1dCjx7QvXtYN1zT7OuPhh+KSJ2tXh2SeOfO\nMGSIZmoWknYIEpGi++ILOP542Gef0BNXEi8sJXIRKaqvvoKTTgr7ZI4cGTaCkMJSIheRolm/Pky5\n37AhzNxs3DjqiMqTZnaKSFFUVcGAAaE2PnasknjUlMhFJG9XXQXvvQfPPBPGi0u0lMhFJC933hkS\n+Msva5x4qVCNXERy9uyzcN55Ye2UffaJOprKoBq5iBTM3Llw9tnw+ONK4qVGMztFJKslS+CEE8J6\n4l27Rh2NZFIiF5EarVkDJ54I558PZ5wRdTRSHdXIRWSLqqrg9NOhSZMw4UezNuufauQiUifXXAMf\nfxyWo1USL11K5CJSrTFjYMQImD5dY8VLnUorIrKZt96C734XnnoqrGgo0cmltKKLnSLyDatWhYWw\nrrtOSTwu1CMXkU3c4dRTYaedwpK0Er0698jNrI2ZTTSzN8xsrpldVk2bhJl9ZmYzU7dBdQ1cRKJx\n442weDHcfnvUkUg+sl3sXAdc4e6zzKwp8JqZ/dPd38xo94K79ylOiCJSH/75T7j1Vpg6VRc346bG\nHrm7L3H3Wan7q4E3gT2qaaqBSSIx9u67Yfr9gw9C69ZRRyP5yvlip5m1BToAUzJecuBIM5ttZuPM\n7KDChScixfbll9C3L/zyl5BIRB2N1EZO48hTZZV/AD9L9czTzQDauPsaM+sFjAHaVXecwYMHb7qf\nSCRI6LdGJHJXXgl77QVXXBF1JAKQTCZJJpN5vSfrqBUzaww8BTzj7rdkPaDZv4BO7r4i43mNWhEp\nMY88Ar/5Dbz2GjRrFnU0Up06T9E3MwOGAfO2lMTNrCWwzN3dzDoTvhxWVNdWRErH22/DJZeENcaV\nxOMtW2mlK3AW8LqZzUw9dzWwJ4C73w30BX5sZuuBNUC/IsUqIgXyxRehLn7ttdCxY9TRSF1pQpBI\nBfrRj0IyHzVKi2GVOq1+KCKbGTECJk+GadOUxMuFeuQiFWTuXOjWDSZOhG9/O+poJBdaNEtENlm9\nOqyj8qc/KYmXG/XIRSqAO5xzDjRuDMOHRx2N5EM1chEBQvKeOTOsoyLlRz1ykTL3+utw7LHw4otw\n4IFRRyP5Uo1cpMKtWhXq4jfdpCReztQjFylT7nDWWdCkCQwbFnU0UluqkYtUsHvvDWWVKZnrlUrZ\nUY9cpAzNng3du8NLL8EBB0QdjdSFauQiFeg//wl18SFDlMQrhXrkImVk4+bJLVrAXXdFHY0Ugmrk\nIhVmyJCwbduoUVFHIvVJPXKRMjFpEpxySri42bZt1NFIoahGLlIhli6Ffv3gvvuUxCuRErlIzK1f\nH5L4+edD795RRyNRUGlFJOYGDoTp02H8eGjYMOpopNB0sVOkzD3xBDzwQNg8WUm8cimRi8TUggVh\ny7axY2GXXaKORqKkGrlIDH32GfTpA3/4AxxxRNTRSNRUIxeJmQ0b4MQTYc894Y47oo5Gik3DD0XK\n0KBB8PnncOutUUcipUI1cpEYeegheOSRsNNP48ZRRyOlosYeuZm1MbOJZvaGmc01s8u20O42M1tg\nZrPNrENxQhWpbNOnw2WXwZgxYS0VkY2y9cjXAVe4+ywzawq8Zmb/dPc3NzYws97Avu6+n5kdDtwJ\ndCleyCKVZ8kSOPlkuPtuaN8+6mik1NTYI3f3Je4+K3V/NfAmsEdGsz7AiFSbKUBzM2tZhFhFKtKX\nX4YkPmBA+CmSKeeLnWbWFugAZO430gpYnPb4A6B1XQMTEaiqgrPPhtat4be/jToaKVU5XexMlVX+\nAfws1TPfrEnG42rHGQ4ePHjT/UQiQSKRyClIkUrkDpdfDsuWhen3DTTGrCIkk0mSyWRe78k6jtzM\nGgNPAc+4+y3VvH4XkHT3h1OP5wNHu/vSjHYaRy6ShxtuCNPvX3wRmjePOhqJSp3HkZuZAcOAedUl\n8ZQngHNS7bsAKzOTuIjkZ8QIuPNOeOYZJXHJrsYeuZkdBbwIvM7X5ZKrgT0B3P3uVLu/AD2Bz4Hz\n3X1GNcdSj1wkB+PHw3nnwcSJcOCBUUcjUculR64p+iIlZNq0sKb42LFw5JFRRyOlQFP0RWLkrbfC\nQljDhimJS36UyEVKwBtvwDHHwPXXh2Qukg+ttSISsdmzoWdP+POf4cwzo45G4kiJXCRCr70Gxx8P\nt98Op54adTQSV0rkIhF59dVQRrnnHvjhD6OORuJMiVwkApMmhXVT7rsv9MhF6kIXO0Xq2XPPwUkn\nwahRSuJSGErkIvXEHYYMgXPPhdGjoUePqCOScqHSikg9WLsWLr4YZs2CV16Btm2jjkjKiXrkIkX2\n8ceQSMCaNTB5spK4FJ4SuUgRTZ0Khx0GJ5wQ9trcbruoI5JypNKKSBG4w9ChYcf7oUPhxBOjjkjK\nmRK5SIG98w5ceCGsWgXJJBx0UNQRSblTaUWkQNavD9Psu3QJwwpfeUVJXOqHeuQiBTB7NlxwATRr\nBlOmwD77RB2RVBL1yEXqYOlS+MUv4Ljj4Mc/hgkTlMSl/imRi9TC++/DpZeGHXzWrAnjwy+4AKzG\n5f9FikOJXCQPb70FAwZAhw6w7bYwbx789a+wxx5RRyaVTDVykSxWrYJx48I48EmT4JJLwsiUHXeM\nOjKRQIlcpBorVsATT4Q1UZJJ6NoVTjkFRo6Epk2jjk7km7T5slS8DRtCD3vmzFDrnjo1bPhw7LEh\neR9/PDRvHnWUUqly2XxZiVzKWlUVfP45fPJJGGGybFm4LV0KH3wAr78ebrvuGurehx4KHTuGtVE0\nnV5KgRK5lLSqKli4EBYv/maC3Xh/9Wr46qvNb+vXh/dv/HXa+HPDhs3bbtgQLkruuuvXt5Ytw8/d\nd4f27eGQQ9TjltJVkERuZsOB44Fl7n5wNa8ngLHAotRTj7n7tdW0UyKvYGvXwty5oXSxsYTx+uvQ\nokVYDTAzye6yC+ywA2y11ea3hg2/HuaX/rNBA9h66y23FYmjQiXy7wKrgZE1JPKfu3ufLMdRIq8w\n6RcMJ06Evff+unzRoYN6wiK5yCWRZx214u4vmVnbbJ+VR1xSxpYsgTFj4LHHwlT17t3htNPCaA8l\nbZHiKMTwQweONLPZwIfAVe4+rwDHlRh54w344x/h6aehd++wG86YMbpgKFIfCpHIZwBt3H2NmfUC\nxgDtqms4ePDgTfcTiQSJRKIAHy9RmjYNrrsurPR3+eVhlmOzZlFHJRJfyWSSZDKZ13tyGrWSKq08\nWV2NvJq2/wI6ufuKjOdVIy8jyWRI4PPnh0WjLrggjA4RkcIqSI08hw9pSRjR4mbWmfDlsCLb+ySe\nliwJU9RnzoTf/AbOOiuMDhGR6GRdNMvMHgJeBvY3s8VmNsDMLjKzi1JN+gJzzGwWcAvQr3jhSlTc\nYfjwMO66XbtQEx8wQElcpBRoQpBktWhR2Lrs009h2LAwfFBE6kcupRUtYytbtGED3HwzdO4M3/9+\nGE6oJC5SerT6oVTr00/hjDPCpgmvvgr77ht1RCKyJeqRy2bmz4fDD4cDDoDnn1cSFyl1SuTyDU8/\nDd/7HgwcCEOGQCP9zSZS8vTfVIAwKuWGG+D222HsWDjiiKgjEpFcKZELa9aECT0LF4ZNFVq1ijoi\nEcmHSisVbtUq6NkzLAH7wgtK4iJxpERewT77DHr0gIMOgvvvh222iToiEakNJfIKtWJFWGL2sMPg\nzjtDj1xE4kn/fSvQ8uVhY+FEAm69VTvoiMSdEnmFWbo0JPDeveHGG5XERcqBEnkF+eijkMRPOw2u\nvVZJXKRcaNGsCvHpp3DUUXDmmXD11VFHIyK5KsjmywUMRok8Il98EUandO4MN90UdTQikg8lcmHD\nBujbNwwtHDVKo1NE4qZedgiS0uUOP/0prF4NjzyiJC5SrpTIy9i114Yp98mkdvIRKWdK5GXq3nvh\nb3+DyZNhhx2ijkZEikk18jL05JNha7YXX4T99os6GhGpC13srEBz5sAxx4R1xTt3jjoaEakr7dlZ\nYZYvhxNPDNPulcRFKod65GVi3bowVvzww+H666OORkQKRaWVCnLJJfDuu2F3n4YNo45GRAqlIKUV\nMxtuZkvNbE4NbW4zswVmNtvMOtQmWKm9oUPDJskPPKAkLlKJcqmR3wf03NKLZtYb2Nfd9wMuBO4s\nUGySg0mTYNCg0BNv1izqaEQkClkTubu/BHxaQ5M+wIhU2ylAczNrWZjwpCbvvx9WMhw5Etq1izoa\nEYlKIUattAIWpz3+AGhdgONKDdauhZNOgiuvhO9/P+poRCRKhZrZmVmIr/aq5uDBgzfdTyQSJBKJ\nAn185bn00tAL//nPo45ERAopmUySTCbzek9Oo1bMrC3wpLsfXM1rdwFJd3849Xg+cLS7L81op1Er\nBXLffWF3n2nToGnTqKMRkWKqrwlBTwDnpD6wC7AyM4lL4cyaBb/8JTz2mJK4iARZSytm9hBwNNDC\nzBYD1wCNAdz9bncfZ2a9zewd4HPg/GIGXMlWrgxri99+Oxx0UNTRiEip0ISgmHAPFzfbtAmJXEQq\ngzaWKCN/+hMsWQKPPhp1JCJSapTIYyCZhJtvDhc3tUGEiGTS6oclbsmSsPP9yJGhrCIikkmJvIRt\n2ABnnAE/+lFY2VBEpDpK5CXs978PP3/722jjEJHSphp5iZowIaxqOGOGVjQUkZqpR16CPv4YzjkH\nRo2C3XaLOhoRKXVK5CVm/fpQF7/oorD3pohINkrkJeb3vw+llEGDoo5EROJCNfIS8txzMGyY6uIi\nkh8l8hLx0Udw7rnw4IPQUttyiEgeVFopAevXQ79+8JOfQLduUUcjInGjRbNKwMCBMHMmjBsHDfTV\nKiJptGhWDDz9dBhmOGOGkriI1I4SeYTeew8GDIDRo2GXXaKORkTiSn3AiHz1FZx2GvziF9C1a9TR\niEicqUYekcsvh0WLYOxYsBqrXyJSyVQjL1GPPRYS+IwZSuIiUnfqkdezBQvgyCPDCJXDDos6GhEp\ndbn0yFUjr0erV8PJJ8PvfqckLiKFox55PXGH/v1hm21g+HCVVEQkN6qRl5AhQ0JZZdIkJXERKays\npRUz62lm881sgZn9qprXE2b2mZnNTN20bl+GiRPhxhvDePFttok6GhEpNzX2yM2sIfAXoDvwITDN\nzJ5w9zczmr7g7n2KFGOsLV4c1hcfNQr22ivqaESkHGXrkXcG3nH3d919HfAwcGI17VQsqMbatXDK\nKXDFFdC9e9TRiEi5ypbIWwGL0x5/kHounQNHmtlsMxtnZgcVMsA4u/RS2HPPMHtTRKRYsl3szGWY\nyQygjbuvMbNewBigXZ0ji7k77oDJk2HKFF3cFJHiypbIPwTapD1uQ+iVb+Luq9LuP2Nmd5jZTu6+\nIvNggwcP3nQ/kUiQSCRqEXLpGz8+bNk2eTJsv33U0YhInCSTSZLJZF7vqXEcuZk1At4CjgU+AqYC\n/dMvdppZS2CZu7uZdQYedfe21RyrIsaRz50bNod4/HE46qiooxGRuKvzOHJ3X29mlwDPAg2BYe7+\nppldlHr9bqAv8GMzWw+sAfoVJPoYWrIETjgBbrlFSVxE6o9mdhbImjWhJ967N1xzTdTRiEi5yKVH\nrkReAFVVcPrpsNVWYby4Lm6KSKFoin49GTQIPv4YJkxQEheR+qdEXkdDh8Kjj8Krr0KTJlFHIyKV\nSIm8Dh56CAYPDmuptGgRdTQiUqmUyGtpzJgw9X7CBGhX8dOfRCRKSuS18OyzcOGF8Mwz8O1vRx2N\niFQ6JfI8vfginHVW6JF36hR1NCIi2uotL1OnQt++8PDD0LVr1NGIiARK5DmaPRt+8IOwTduxx0Yd\njYjI15TIc/DSS9CjB/zlL2EKvohIKVEiz2L06LA5xKhRcOqpUUcjIrI5XeyswR13wLXXhmVpO3aM\nOhoRkeopkVfDPUy7//vfw673e+8ddUQiIlumRJ5h3Tq46KKwrvjkybDLLlFHJCJSMyXyNMuWwdln\nQ6NGYdr9dttFHZGISHa62Jny/PPQoUOohY8ZoyQuIvFR8T3ydevCRhAjRoRb9+5RRyQikp+KTuTv\nvgv9+8OOO8LMmbDrrlFHJCKSv4osrbjDI49A585hyv1TTymJi0h8VVyPfMYMuPLKcGFz3Dj4znei\njkhEpG4qpkf+4Ydw3nlhc+TTTw9rpyiJi0g5KPtEvnp1uJjZvj3svju8/TZcfHEYYigiUg7KNp29\n9x4MGwb33guJRCip7LVX1FGJiBRe1h65mfU0s/lmtsDMfrWFNrelXp9tZh0KH2Zu1q2Dxx8P5ZOO\nHWHlSnjuOXjwQSVxESlfNSZyM2sI/AXoCRwE9DezAzPa9Ab2dff9gAuBO4sUa7XWrYNXXoGrrw7J\n+uabw5DCDz6A224r3FZsyWSyMAeKQJxjB8UfNcVf+rL1yDsD77j7u+6+DngYODGjTR9gBIC7TwGa\nm1nLgkeasn49TJkC118PPXvCzjvDT34SEvqECWHt8LPPhm22KeznxvmXIc6xg+KPmuIvfdlq5K2A\nxWmPPwAOz6FNa2BpbYP68sswPHDhQli0KNw23p8/H9q2DXXviy8OZZOddqrtJ4mIxF+2RO45Hsdy\neV+vXlBV9c3b2rWwalW4rV4dfgK0aAH77BOWkN1nHzj++HB///1DL1xERAJz33KuNrMuwGB375l6\nPBCocvcb0trcBSTd/eHU4/nA0e6+NONYuX4piIhIGnfP7Cx/Q7Ye+XRgPzNrC3wEnA70z2jzBHAJ\n8HAq8a/MTOK5BCIiIrVTYyJ39/VmdgnwLNAQGObub5rZRanX73b3cWbW28zeAT4Hzi961CIiskmN\npRURESl9RZ+in8uEolJlZsPNbKmZzYk6ltowszZmNtHM3jCzuWZ2WdQx5cPMmpjZFDObZWbzzOyP\nUceULzNraGYzzezJqGOpDTN718xeT/0bpkYdTz7MrLmZ/cPM3kz9/nSJOqZcmdn+qXO+8fZZTf9/\ni9ojT00oegvoDnwITAP6u/ubRfvQAjKz7wKrgZHufnDU8eTLzHYDdnP3WWbWFHgN+GFczj+AmW3r\n7mvMrBEwCbjK3SdFHVeuzOznQCdge3fvE3U8+TKzfwGd3H1F1LHky8xGAC+4+/DU78927v5Z1HHl\ny8waEPJWE/2+AAACU0lEQVRnZ3dfXF2bYvfIc5lQVLLc/SXg06jjqC13X+Lus1L3VwNvAntEG1V+\n3H1N6u5WhOs0sUkoZtYa6A3cy+ZDdOMkdrGbWTPgu+4+HML1vjgm8ZTuwMItJXEofiKvbrJQqyJ/\nplQjNfKoAzAl2kjyY2YNzGwWYYLZRHefF3VMeRgC/AKoijqQOnBggplNN7P/jjqYPHwL+MTM7jOz\nGWY21My2jTqoWuoHPFhTg2Incl1JLQGpsso/gJ+leuax4e5V7n4oYbbw98wsEXFIOTGzE4Bl7j6T\nGPZo03R19w5AL+CnqXJjHDQCOgJ3uHtHwoi6X0cbUv7MbCvgB8Dfa2pX7ET+IdAm7XEbQq9c6omZ\nNQYeA0a5+5io46mt1J/FTwNx2Q7kSKBPqsb8EHCMmY2MOKa8ufvHqZ+fAI8TyqVx8AHwgbtPSz3+\nByGxx00v4LXU+d+iYifyTROKUt8spxMmEEk9MDMDhgHz3P2WqOPJl5m1MLPmqfvbAMcBM6ONKjfu\nfrW7t3H3bxH+NP4/dz8n6rjyYWbbmtn2qfvbAT2AWIzgcvclwGIza5d6qjvwRoQh1VZ/QkegRkXd\nWGJLE4qK+ZmFZGYPAUcDO5vZYuC37n5fxGHloytwFvC6mW1MgAPdfXyEMeVjd2BE6qp9A+B+d38+\n4phqK45lxpbA46E/QCPgAXd/LtqQ8nIp8ECqE7mQmE1WTH15dgeyXpvQhCARkZgr+z07RUTKnRK5\niEjMKZGLiMScErmISMwpkYuIxJwSuYhIzCmRi4jEnBK5iEjM/T8bTzxBOTN6TgAAAABJRU5ErkJg\ngg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "x = np.linspace(0, 2 * np.pi)\n", "p = plt.plot(x, Y(x))\n", "t = plt.title(r'$Y(x) = \\int_0^x sin^2(\\theta) d\\theta$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 数值积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "数值积分:\n", "\n", "$$F(x) = \\lim_{n \\rightarrow \\infty} \\sum_{i=0}^{n-1} f(x_i)(x_{i+1}-x_i) \n", "\\Rightarrow F(x) = \\int_{x_0}^{x_n} f(x) dx$$\n", "\n", "导入贝塞尔函数:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.special import jv" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def f(x):\n", " return jv(2.5, x)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczXX///HHy1K2EqlMsu+T0W9sKWFs0UxX8lOWyKip\nK1kGaZkuhK6ISiHL2AalqMs6XBVizkiLRk0MMxiFLFGSNNbB+/vHDM3FDDPnzJz3WV73283t9vmc\n8zmf99O51cv7vD+fz/stxhiUUkr5h0K2AyillHIfLfpKKeVHtOgrpZQf0aKvlFJ+RIu+Ukr5ES36\nSinlR1wu+iLSQUS2i0iqiLyUzfshIvKniCRm/hnmaptKKaWcU8SVD4tIYWAy0BY4ACSISKwxJuWy\nQ+ONMQ+50pZSSinXudrTbwLsMsbsMcakAwuBjtkcJy62o5RSKh+4WvQrAPuy7O/PfC0rA9wrIptF\n5BMRCXSxTaWUUk5yaXiHjIJ+Ld8DFY0xJ0XkAWAZUMvFdpVSSjnB1aJ/AKiYZb8iGb39S4wxf2XZ\n/lREpopIWWPM0azHiYhOAqSUUk4wxuR6CN3V4Z1NQE0RqSIi1wFdgdisB4jIbSIimdtNALm84F9k\njNE/xjBixAjrGTzlj34X+l3od3H1P3nlUk/fGHNORPoDq4DCwGxjTIqIPJP5/nTgEeBZETkHnAS6\nudKmUkop57k6vIMx5lPg08tem55lewowxdV2lFJKuU6fyPVAISEhtiN4DP0u/qbfxd/0u3CeODMm\nVBBExHhKFqWU8hYignHjhVyllFJeRIu+Ukr5ES36SinlR1y+e0ep3Dhx4gQHDhzg4MGDVKlShSpV\nqtiOpJRf0qKv8tWZM2f48MMP2bBhAwcOHGD//v0cOHCAU6dOcccddxAQEEBqaiolSpSgdevWtG7d\nmlatWhEQEGA7ulJ+Qe/eUfni6NGjREdH8+6773LXXXfRqVMnKlasSIUKFbjjjjsoW7YsmQ9mY4wh\nJSWFdevWsW7dOhwOBwEBAbRp04bnnntOfwUolQd5vXtHi75yye7du5kwYQLvv/8+Dz30EEOGDCEo\nKChP5zh//jybN29myZIlREdH869//YvIyEiKFNEfokpdixZ95Ra7du1i6NChrF27lqeeeooBAwZQ\nocLls2rnXWpqKn369OHYsWPMmDGDhg0b5kNapXyX3qevCtySJUu49957CQ4OZvfu3YwdOzZfCj5A\nzZo1+fzzz4mMjCQ0NJQhQ4aQlpaWL+dWSmnRV3mQnp7O888/z3PPPcd///tfoqKiuOGGG/K9HREh\nPDycrVu38ttvv1GvXj0++eSTfG9HKX+kwzsqVw4ePEjXrl254YYbeP/997n55pvd1vaaNWt48skn\nGTZsGM8884zb2lXKG+jwjsp369ato1GjRnTo0IGVK1e6teADtGvXDofDwejRo4mOjnZr20r5Gr09\nQuXowoULjB07lnfffZf58+fTpk0ba1mqV69OXFwcrVu3BqBPnz7WsijlzbToq2wZY4iMjCQhIYFN\nmzbl24VaV1SvXp1169bRunVrjDE8++yztiMp5XVcHt4RkQ4isl1EUkXkpasc11hEzonI/3e1TVWw\njDG8+OKLbNy4kdWrV3tEwb/oYuEfO3Ys06ZNsx1HKa/jUk9fRAoDk4G2ZCySniAiscaYlGyOGwd8\nBuT6goOyY9SoUaxevZq4uDhKly5tO84VLg71tGrVCmMMffv2tR1JKa/h6vBOE2CXMWYPgIgsBDoC\nKZcdNwBYBDR2sT1VwMaNG8dHH31EfHw8ZcuWtR0nR9WqVbtU+EVEh3qUyiVXi34FYF+W/f3A3VkP\nEJEKZPxD0JqMoq/3ZXqoSZMmMWPGDNavX8+tt95qO841XSz89957L/Xq1aN58+a2Iynl8Vwt+rkp\n4BOAKGOMkYwZt3Ic3hk5cuSl7ZCQEF0H041mzpzJ+PHjiY+P96gx/GupVq0aMTExPPbYYyQmJlKu\nXDnbkZQqUA6HA4fD4fTnXXo4S0SaAiONMR0y918GLhhjxmU55if+LvTlgJPA08aY2MvOpQ9nWTJ/\n/nyioqKIi4ujZs2atuM45cUXXyQ5OZnY2FgKFdLHT5T/cOuEayJSBNgBtAEOAt8C3S+/kJvl+DnA\nCmPMkmze06JvQVxcHN27d2fdunUEBgbajuO09PR0WrRoQefOnXn++edtx1HKbfJa9F0a3jHGnBOR\n/sAqoDAw2xiTIiLPZL4/3ZXzq4L1888/89hjjzF//nyvLvgARYsWZeHChTRp0oT77ruPpk2b2o6k\nlEfSuXf81OnTp2nevDldunThhRdesB0n3yxbtoxBgwaRmJhImTJlbMdRqsDpfPrqmowxREREcOLE\nCRYuXHhpRStfMXDgQH7++WeWLFnic383pS6nE66pa4qOjiYhIYGYmBifLIpvvPEG+/btY/Lkybaj\nKOVxtKfvZ7766is6derEl19+SY0aNWzHKTA//vgjTZs25bPPPtPVt5RP056+ytEvv/xCly5dmDNn\njk8XfMiYqmHSpEn06tWL9PR023GU8hha9P3E2bNneeSRR3jmmWcIDQ21HcctunXrRsWKFZk0aZLt\nKEp5DB3e8RP9+/dn3759LF261K8eXkpNTeWee+5h8+bNXvWksVK5pcM76gqxsbF88sknvPfee35V\n8CFjofU+ffowZMgQ21GU8gja0/dxhw4dIjg4mEWLFtGsWTPbcaw4efIkd955J7NmzbK6+pdSBUF7\n+uqSi/fjR0RE+G3BByhRogQTJkygf//+nD171nYcpazSou/DoqOj+fXXXxkxYoTtKNY99NBDVKtW\njXfeecd2FKWs0uEdH7Vjxw7uu+8+NmzYQO3atW3H8Qg//vgjd999N4mJiVSsWNF2HKXyhQ7vKNLT\n0+nRowevvvqqFvwsqlevTv/+/Rk8eLDtKEpZoz19HzRs2DASExNZuXKlT06z4IpTp05Rr149pk6d\nSvv27W3HUcplOuGan9uwYQOPPvooP/zwA7fddpvtOB7pv//9L4MHDyYpKYnrr7/edhylXKLDO37s\n+PHj9OrVi+nTp2vBv4qwsDDq1q3LhAkTbEdRyu20p+9DIiIiKFy4MDNmzLAdxeNt376d5s2bk5qa\nyk033WQ7jlJOc3tPX0Q6iMh2EUkVkZeyeb+jiGwWkUQR+U5EWrvaprrSqlWrWLt2LePHj7cdxSvU\nqVOHsLAw3n77bdtRlHIrV9fILUzGGrltgQNAApetkSsiJY0xJzK3g4ClxpgrpnjUnr7zjh8/TlBQ\nEDNnzuT++++3Hcdr7N69m0aNGrFjxw7KlStnO45STnF3T78JsMsYs8cYkw4sBDpmPeBiwc9UCjji\nYpvqMlFRUbRt21YLfh5VrVqVrl27Mm7cONtRlHIblxZGByoA+7Ls7wfuvvwgEXkYeB0IALQy5aP4\n+HhiY2PZunWr7SheaejQoQQFBTF48GBuv/1223GUKnCuFv1cjccYY5YBy0SkOfA+kO0TQyNHjry0\nHRISQkhIiIvxfNvJkyeJiIhg6tSpejHSSRUqVOCJJ55gzJgxuryi8goOhwOHw+H0510d028KjDTG\ndMjcfxm4YIzJ8feyiPwINDHG/H7Z6zqmn0fPP/88Bw4cYMGCBbajeLXffvuNOnXq8N1331GlShXb\ncZTKE7c+nCUiRci4kNsGOAh8y5UXcqsDPxljjIg0AP5jjKmezbm06OfBxo0b6dixI0lJSdxyyy22\n43i9YcOGcfDgQWJiYmxHUSpP8lr0XRreMcacE5H+wCqgMDDbGJMiIs9kvj8d6Az0EpF0IA3o5kqb\nCs6cOcOTTz7JxIkTteDnk+eff56aNWuyc+dOatWqZTuOUgVGH87yQsOHDycpKYmlS5fq3Dr5aMyY\nMSQlJelwmfIqOveOj9u8eTPt2rXjhx9+0LtN8llaWho1atRg9erV1K9f33YcpXJF597xYefOnSMi\nIoKxY8dqwS8ApUqV4qWXXuKVV16xHUWpAqNF34tMmjSJ0qVL88QTT9iO4rOeffZZNm3aREJCgu0o\nShUIHd7xErt376Zx48Z888031KhxxSwWKh9NmjQJh8PBkiVLbEdR6pp0TN8HGWN44IEHCAkJISoq\nynYcn3fixAmqVq3KF198oSuPKY+nY/o+6MMPP+TQoUMMGTLEdhS/ULJkSfr168ebb75pO4pS+U57\n+h7uyJEj1KtXjxUrVtC4cWPbcfzG77//Ts2aNUlKSqJChQq24yiVIx3e8THh4eGULVuWd955x3YU\nvzNo0CCKFi2qPX7l0bTo+5A1a9bw1FNPsW3bNkqVKmU7jt/5+eefCQ4OZteuXZQpU8Z2HKWypWP6\nPuLkyZP06dOHadOmacG3pFKlSoSFhTFt2jTbUZTKN9rT91Avvvgi+/bt0ykBLNu6dStt27Zl9+7d\nFC9e3HYcpa6gwzs+IDExkfbt25OUlMRtt91mO47f+8c//kFoaCjPPvus7ShKXUGLvpc7d+4cTZo0\nITIykt69e9uOo4ANGzYQHh7Ojh07KFLE1XWHlMpfOqbv5caPH0+5cuUIDw+3HUVluu+++wgICGDx\n4sW2oyjlMu3pe5DU1FTuueceEhISqFq1qu04KosVK1bwyiuv8P333+t01sqjaE/fS124cIGnn36a\noUOHasH3QGFhYZw9e5Y1a9bYjqKUS1wu+iLSQUS2i0iqiLyUzfs9RGSziGwRkS9FRCcqz8asWbM4\ndeoUkZGRtqOobBQqVIiXXnqJceNyXP5ZKa/g6hq5hclYI7ctcABI4Mo1cu8Bko0xf4pIBzIWUm+a\nzbn8dnjn4MGD3HXXXaxbt46goCDbcVQO0tPTqV69OkuXLqVhw4a24ygFuH94pwmwyxizxxiTDiwE\nOmY9wBjztTHmz8zdjcAdLrbpU4wx9OvXj2effVYLvocrWrQo/fv3Z+LEibajKOU0V+8/qwDsy7K/\nH7j7KsdHAJ+42KZPWbx4MTt27GDhwoW2o6hceOqpp6hevTq//PILAQEBtuMolWeuFv1cj8eISCvg\nSaBZTseMHDny0nZISAghISEuRPN8R48eJTIykkWLFnH99dfbjqNyoWzZsnTr1o3o6GhGjRplO47y\nQw6HA4fDwYULF1i5cmWeP+/qmH5TMsboO2TuvwxcMMaMu+y4+sASoIMxZlcO5/K7Mf0nn3ySkiVL\n8u6779qOovIgJSWFkJAQ9u7dS7FixWzHUX5q8ODBpKSksGrVKreO6W8CaopIFRG5DugKxGY9QEQq\nkVHwe+ZU8P3R6tWrWbt2LWPGjLEdReVR3bp1CQ4O1iE5Zc3777/PihUrnJqby+WHs0TkAWACUBiY\nbYx5XUSeATDGTBeRWUAn4OfMj6QbY5pkcx6/6ekfO3aMoKAg5syZQ9u2bW3HUU747LPPiIqKIjEx\nUR/WUm713Xff0aFDB+Li4qhXr57OveMNevfuTcmSJZkyZYrtKMpJFy5cIDAwkOnTp9OyZUvbcZSf\n+PXXX2ncuDFvv/02nTt3BvJ+y6bOHuVmsbGxbNiwgR9++MF2FOWCQoUKMXDgQCZMmKBFX7lFeno6\njz76KI8//vilgu8M7em70e+//05QUBAfffQRzZs3tx1HuejEiRNUrlyZb7/9lmrVqtmOo3xcZGQk\nP/74I7GxsRQuXPjS6zr3jgfr168f3bp104LvI0qWLMmTTz7J5MmTbUdRPm7u3Ll89tlnfPDBB/9T\n8J2hPX03+fjjj3nllVdITEzUFZh8yMV1dPfs2cMNN9xgO47yQQkJCYSGhhIfH09gYOAV72tP3wMd\nPnyYyMhI5s2bpwXfx1SqVInWrVszd+5c21GUD0pLS6Nr165Mnz4924LvDO3pFzBjDJ06dSIwMFDv\nyfdRX375Jb1792bHjh0UKqT9KJV/+vXrx8mTJ5kzZ06Ox+jdOx5m/vz5/PTTT3z00Ue2o6gCcu+9\n91K6dGk++eQTHnzwQdtxlI+Ii4tj+fLlJCUl5et5tVtSgPbu3cuQIUOYN2+ezq3jw0SEQYMGMWHC\nBNtRlI9IS0sjIiKC6dOnU6ZMmXw9tw7vFJD09HRatGhB586def75523HUQXs7NmzVK5cmbVr1+bb\n2KvyXwMGDOCvv/7K1bUifSLXQ0RFRbFlyxZWrlyp47x+4pVXXuHo0aN6C6dyicPhoGfPniQlJeWq\nl69F3wOsWrWKiIgIEhMTueWWW2zHUW5y4MABgoKC2LNnDzfeeKPtOMoLnThxgvr16zNhwgT+8Y9/\n5OozesumZQcPHqR379588MEHWvD9TIUKFWjbti3vvfee7SjKS7388ss0a9Ys1wXfGdrTz0fnz5+n\nXbt2tGzZkhEjRtiOoyyIj4+nT58+JCcn6+ybKk/i4+N57LHH2Lp1a54u3mpP36IxY8ZgjGHYsGG2\noyhLWrRoQZEiRVi3bp3tKMqLnDhxgoiICKKjo/P9bp3LaU8/n6xfv56uXbvy3Xffcfvtt9uOoyyK\njo5m1apVLF261HYU5SWGDBnC4cOHmT9/fp4/6/YLuSLSgb8XUZmVzVKJdYA5QDAw1BgzPofzeG3R\nP3LkCMHBwcycOZMOHTrYjqMsS0tLo3LlyiQmJlKpUiXbcZSHS05OpmXLlmzbto1bb701z5936/CO\niBQGJgMdgECgu4jUveyw34EBwFuutOWpzp8/T3h4ON27d9eCrwAoVaoUjz/+ONHR0bajKA9njGHg\nwIEMGzbMqYLvDFfH9JsAu4wxe4wx6cBCoGPWA4wxvxljNgHpLrblkaKiojh16hSjR4+2HUV5kL59\n+zJr1ixOnz5tO4ryYMuXL+fgwYP07dvXbW26WvQrAPuy7O/PfM0vzJo1i+XLl7No0SKKFi1qO47y\nILVq1SI4OJj//Oc/tqMoD3X69Gmee+45Jk6c6Nb64WrR985B+HwQFxfH0KFDWblyJWXLlrUdR3mg\nfv366dO5Kkfjx48nODiYtm3burVdV2fZPABUzLJfkYzevlNGjhx5aTskJISQkBBnT1Wgdu7cSbdu\n3ViwYAG1atWyHUd5qLCwMCIjI0lISKBx48a24ygPsm/fPt5++202bdqU5886HA4cDofTbbt0946I\nFAF2AG2Ag8C3QHdjTEo2x44E/vL2u3eOHj1K06ZNefHFF3nqqadsx1Ee7o033iA5OVkXWVH/o3v3\n7tSsWZNXX33V5XPZuGXzAf6+ZXO2MeZ1EXkGwBgzXUTKAwnAjcAF4C8g0BiTdtl5PL7onz17lg4d\nOtCgQQPeessnb0ZS+ezIkSPUqFGD1NRUnZZDARnP9PTs2ZPt27dTokQJl8+nE64VEGMM//znPzl0\n6BDLli1zeXFi5T+eeOIJateuTVRUlO0oyrJz587RsGFDhg4dSpcuXfLlnDoNQwF56623+Pbbb/nw\nww+14Ks86devH9OmTeP8+fO2oyjLZs6cSZkyZXj00UetZdCinwvvvPMO06ZNY+XKldxwww224ygv\n06hRIwICAli5cqXtKMqi33//nREjRjBp0iSrk/Fp0b+Gt956iylTpuBwOKhYseK1P6BUNvr168eU\nKVNsx1AW/fvf/+aRRx6hfv36VnPomP5VvPHGG8ycOZO4uDjuuOMO23GUFzt9+jSVKlViw4YNepuv\nH9q9ezeNGjUiOTmZ2267LV/PrWP6+eT1119n1qxZOBwOLfjKZcWKFSMiIoKpU6fajqIsGD58OJGR\nkfle8J2hPf1sjB49mvfee4+4uDidJlnlm71799KgQQP27t1LqVKlbMdRbpKYmEhoaCg7d+4skGuC\n2tN30auvvsr8+fNxOBxa8FW+qly5Ms2bN+eDDz6wHUW5UVRUFMOHD/eYm0C06Gc6d+4cL7zwAgsX\nLiQuLo6AgADbkZQPunhB11N+1aqC9fnnn/PTTz/x9NNP245yiRZ94PDhw9x///1s3ryZ9evXU758\neduRlI9q06YNZ86cYcOGDbajqAJ24cIFXnrpJcaMGeNRs/D6fdH/6quvaNSoEc2aNePTTz+lXLly\ntiMpH1aoUCG9fdNPfPTRRxQuXJhHHnnEdpT/4bcXco0xTJ48mddee42YmBjCwsLc1rbyb3/++SdV\nqlQhOTlZhxF91NmzZ6lTpw4xMTEFPluwXsjNhbS0NHr06EFMTAxff/21FnzlVqVLl6Zr167MnDnT\ndhRVQKZPn06dOnU8cnp4v+vpJyYm0rNnT+6++26mTJlC8eLFC7xNpS6XlJREhw4d2LNnj0eN9yrX\nHT9+nFq1arF69Wq3PH2rPf0c7N+/n969exMaGsqLL75ITEyMFnxlTVBQEDVq1GDZsmW2o6h89tZb\nb9G+fXvr0y3kxOeL/l9//cXw4cO56667qFChAjt27CA8PNx2LKX0gq4POnToEFOmTMmXxVEKis8W\n/XPnzjFjxgxq167N3r17SUxMZPTo0dx44422oykFQKdOndi5cydbt261HUXlk9GjRxMeHk7lypVt\nR8lRfqyc1YG/V86aZYwZl80xk4AHgJNAb2NMYjbH5MuY/m+//cZ//vMfpk6dSrly5Rg/fjwNGzZ0\n+bxKFYRRo0Zx6NAhpk2bZjuKctHFaTZSUlK49dZb3dauW1fOEpHCZKyR25aMRdITuGyNXBEJBfob\nY0JF5G5gojGmaTbncrroHz9+nKVLl7JgwQK++eYbQkNDCQ8P5/7777c6b7VS1/LLL79w55138tNP\nP3HTTTfZjqNcEBERQUBAAK+99ppb281r0S/iYntNgF3GmD2ZjS8EOgJZF0Z/CJgHYIzZKCI3icht\nxpjDzjZ64cIF9uzZw6ZNm/j4449Zs2YNISEh9O7dm8WLF1OyZEnn/0ZKuVFAQAChoaHMnj2bIUOG\n2I6jnLRz505iY2NJTU21HeWaXC36FYB9Wfb3A3fn4pg7gGyLvjGG06dPc+rUKU6ePMmJEydITU1l\n27ZtbNu2jeTkZFJSUihbtiz169enU6dOl5YgU8obDRw4kC5dujBo0CBditNLjRgxgsGDB3vFrzVX\ni35ux2Mu/+mR7eeKFy/OmTNnuP766ylevDjFixenRIkSVK9encDAQFq2bEnfvn0JDAzUC7LKZzRu\n3Jjy5cuzcuVKOnbsaDuOyqMtW7YQFxfnNQ/buVr0DwBZ1xCsSEZP/mrH3JH52hWee+45ihQpgogQ\nEhLikU+zKVUQIiMjmThxohZ9LzR8+HCioqLctkaCw+HA4XA4/XlXL+QWIeNCbhvgIPAtV7+Q2xSY\nkN8XcpXydmfPnqVq1ap89tlnBAUF2Y6jcmnjxo088sgjpKamUqxYMSsZ3PpErjHmHNAfWAUkAx8Z\nY1JE5BkReSbzmE+An0RkFzAd6OtKm0r5ouuuu45nn32Wd99913YUlQfDhg1j+PDh1gq+M/xu7h2l\nPNWvv/5K7dq12bVrFzfffLPtOOoaHA4HTz31FCkpKVbnT9K5d5TyUrfeeisdO3Zk1qxZtqOoazDG\nMHToUEaOHOl1E+Zp0VfKg0RGRjJlyhTOnTtnO4q6ik8//ZRjx47RvXt321HyTIu+Uh6kQYMGVK5c\nmeXLl9uOonJw4cIFhg0bxr///W+vfK5Ci75SHubi7ZvKMy1ZsoRChQrRqVMn21GcohdylfIw6enp\nVKtWjdjYWIKDg23HUVmcP3+eoKAg3n77bTp06GA7DqAXcpXyekWLFqVv3756+6YH+uCDD7j55ptp\n37697ShO056+Uh7oyJEj1KxZk507d3LLLbfYjqP4e7HzuXPn0qJFC9txLtGevlI+oFy5cnTu3Jno\n6GjbUVSmmJgYatas6VEF3xna01fKQyUnJ9O6dWv27NnjVU98+qJTp05Rs2ZNli5dSuPGjW3H+R/a\n01fKRwQGBtKoUSPee+8921H83rRp02jcuLHHFXxnaE9fKQ8WHx/P008/TUpKilfeE+4L/vrrL2rW\nrMnnn39OvXr1bMe5gvb0lfIhLVq0oEyZMsTGxtqO4rcmTpxImzZtPLLgO0N7+kp5uEWLFvH222/z\n1Vdf2Y7id/744w9q1arF119/TY0aNWzHyZb29JXyMZ06deLXX3/lyy+/tB3F77z55ps8/PDDHlvw\nnaE9faW8wNSpU1m9ejXLli2zHcVvHD58mMDAQBITE6lUqZLtODnKa09fi75SXuDkyZNUrVqV+Ph4\n6tSpYzuOXxg0aBDGGI+fB8ltRV9EygIfAZWBPUAXY8yxbI6LAcKAX40xOa4Dp0VfqasbNWoU+/fv\n95oFuL3Z3r17adCgAdu2baN8+fK241yVO4v+G8ARY8wbIvISUMYYE5XNcc2BNOA9LfpKOe/IkSPU\nqlWL5ORkjy9E3i48PJzKlSvz6quv2o5yTe4s+tuBlsaYwyJSHnAYY7L93SkiVYAVWvSVck2/fv24\n6aabGD16tO0oPmvLli20a9eO1NRUbrzxRttxrsmdRf8PY0yZzG0Bjl7cz+bYKmjRV8plP/74I02b\nNmX37t2UKlXKdhyfFBYWRvv27YmMjLQdJVfyWvSLXONka4DsfkcOzbpjjDEi4nLFHjly5KXtkJAQ\nQkJCXD2lUj6levXqtGrVitmzZzNw4EDbcXyOw+EgJSWFJUuW2I6SI4fDgcPhcPrzrg7vhBhjDolI\nABCnwztKFbyEhAQeffRRUlNTvW5Rbk9mjKFp06YMHDiQxx57zHacXHPnw1mxQHjmdjigNxAr5QaN\nGzematWqLFy40HYUn7J48WLOnj1Lt27dbEcpUK7esvkxUIkst2yKyO3ATGNMWOZxC4CWwM3Ar8Ar\nxpg52ZxPe/pK5VJcXBz//Oc/SUlJoUiRq47SqlxIT0/nzjvvZPLkydx///224+SJPpyllJ8ICQnh\niSeeIDw8/NoHq6uKjo5m0aJFrFmzhoz7UryHFn2l/ER8fDwRERGkpKTo2L4L0tLSqFWrFrGxsTRq\n1Mh2nDzTCdeU8hMtW7akcuXKvP/++7ajeLUJEybQokULryz4ztCevlJebMOGDfTq1YsdO3Zob98J\nv/32G3Xq1GHjxo1eO5Om9vSV8iP33XcfNWrUYO7cubajeKXXXnuN7t27e23Bd4b29JXycl9//TXd\nu3dn587fphsOAAANG0lEQVSdXHfddbbjeI3k5GRatGhBcnIyt956q+04TtOevlJ+5p577qFu3brE\nxMTYjuI1jDEMGjSIYcOGeXXBd4b29JXyAd9++y2PPPIIqampXH/99bbjeLzly5fz8ssvs3nzZq+/\nFqI9faX8UJMmTQgKCmLWrFm2o3i806dPM3jwYCZOnOj1Bd8Z2tNXykds2rSJhx9+mF27dlGsWDHb\ncTzWmDFjSEhIYOnSpbaj5At9OEspP/bQQw/Rrl07BgwYYDuKR9q/fz933XUXCQkJVKtWzXacfKFF\nXyk/lpiYSFhYGLt27aJEiRK243icHj16ULVqVV577TXbUfKNFn2l/FzXrl0JDAxkxIgRtqN4lA0b\nNtC9e3e2b99OyZIlbcfJN1r0lfJzFxf1/v7776lcubLtOB7h/PnzNG7cmBdeeIHu3bvbjpOv9O4d\npfxc5cqViYyM5IUXXrAdxWPMnj2bkiVL+vxc+bmhPX2lfNCpU6eoW7cuc+fO9ftlR//44w/q1q3L\np59+SnBwsO04+U6Hd5RSACxatIhXX32V77//3q8XWunbty8XLlwgOjradpQC4dbhHREpKyJrRGSn\niKwWkZuyOaaiiMSJyDYR2Soi3rHEvFJernPnzpQrV44ZM2bYjmJNXFwcK1asYOzYsbajeAyXevoi\n8gZwxBjzhoi8BJQxxkRddkx5oLwx5gcRKQV8BzxsjEm57Djt6SuVz5KSkmjTpg0pKSncfPPNtuO4\nVVpaGvXr1+fdd98lLCzMdpwC49bhHRHZDrQ0xhzOLO4OY0yda3xmGfCuMWbtZa9r0VeqAAwYMIDz\n588zdepU21Hcqn///qSlpfn8tNPuLvp/GGPKZG4LcPTifg7HVwHigTuNMWmXvadFX6kCcPToUerW\nrcvq1au56667bMdxC4fDQc+ePUlKSqJMmRxLkk/Ia9G/5tUdEVkDlM/mraFZd4wxRkRyrNqZQzuL\ngIGXF/yLRo4ceWk7JCTE7+86UCo/lC1bllGjRhEZGYnD4fC6hb/zKi0tjYiICKZPn+6TBd/hcOBw\nOJz+fH4M74QYYw6JSAAQl93wjogUBVYCnxpjJuRwLu3pK1VAzp8/T8OGDfnXv/5Fly5dbMcpUAMG\nDOD48ePMmzfPdhS3cPfwzhvA78aYcSISBdyUzYVcAeZlHjf4KufSoq9UAfriiy/o0aMHSUlJlC5d\n2nacAhEfH3/p7+iLvfzsuLvolwU+BioBe4AuxphjInI7MNMYEyYi9wHrgS3AxcZeNsZ8dtm5tOgr\nVcD69u3LsWPH+OCDD3xumOfEiRPUr1+fiRMn8uCDD9qO4zb6cJZSKkcnT56kcePGREVF8fjjj9uO\nk68iIyP5888//WZY5yIt+kqpq9qyZQtt2rThm2++oXr16rbj5It169bRq1cvvxrWuUgnXFNKXVX9\n+vUZNmwYjz32GOnp6bbjuGzPnj306NGDuXPn+l3Bd4b29JXyQ8YYwsLCCA4OZvTo0bbjOO3EiRM0\na9aM3r17M2jQINtxrNDhHaVUrhw+fJjg4GA+/PBDr3wmxhhD165dKVGiBHPmzPG5C9O5pcM7Sqlc\nue2224iJiaFXr14cPXrUdpw8e/3119m7dy/R0dF+W/CdoT19pfzcoEGD2LdvH4sWLfKa4rly5Uqe\neeYZEhISuP32223HsUp7+kqpPBk7diy7du1i1qxZtqPkSkpKCk8++SSLFi3y+4LvDP9dWUEpBUCx\nYsVYsGABrVq1okKFCoSGhtqOlKNjx47RsWNHxo0bxz333GM7jlfS4R2lFADffPMNDz30EAsWLKBN\nmza241zh/PnzPPjgg9SqVYuJEyfajuMx9O4dpZTTvvjiCzp37szixYtp3ry57TiXnD59mp49e5KW\nlsaKFSsoWrSo7UgeQ8f0lVJOa968OQsWLKBz58588803tuMAcPz4cR544AEKFSrE8uXLteC7SIu+\nUup/tGnThnnz5tGxY0e+//57q1kOHz5MSEgIgYGBLFiwgOuvv95qHl+gRV8pdYUHHniAGTNmEBoa\nSlJSkpUMP/30E82aNaNjx45MnjyZwoULW8nha/TuHaVUtjp27MiZM2do3749a9eupW7dum5r+4cf\nfiAsLIzhw4fTp08ft7XrD7Snr5TKUZcuXXjzzTdp3rw5kyZN4vz58wXepsPh4P7772fixIla8AuA\n3r2jlLqmHTt28PTTT5Oens7s2bMJDAzM9zaOHj3KmDFjmDdvHh9//DGtWrXK9zZ8kdvu3hGRsiKy\nRkR2ishqEbkpm2OKichGEflBRJJF5HVn21NK2VO7dm0cDgfh4eG0bNmSUaNGcfbs2Xw595kzZxg/\nfjy1a9fmr7/+YsuWLVrwC5ArwztRwBpjTC1gbeb+/zDGnAZaGWP+H1AfaJW5fKJSyssUKlSIPn36\nkJiYyHfffUeDBg1cuq3zwoULfPjhh9SpU4f4+HjWr1/P9OnTCQgIyMfU6nJOD++IyHagpTHmsIiU\nBxzGmDpXOb4EEA+EG2OSs3lfh3eU8hLGGD7++GMGDRpEw4YNCQkJoUWLFgQHB1/zPvojR46wceNG\nRowYQaFChXjzzTdp2bKlm5L7Hrc9kSsifxhjymRuC3D04v5lxxUCvgeqA9OMMS/mcD4t+kp5mWPH\njrFmzRrWr1/P+vXr2b17N02bNqVFixbcd999nDlzhpSUlP/5k56ezp133smAAQPo0qULhQrp/SSu\nyNeiLyJrgPLZvDUUmJe1yIvIUWNM2aucqzSwCogyxjiyed+MGDHi0n5ISIhXLuyglD87evQoX375\nJevXr2fDhg2UKFGCunXr/s+f8uXLe80Uzp7I4XDgcDgu7Y8aNcptPf3tQIgx5pCIBABxVxveyfzM\ncOCUMeatbN7Tnr5SSuWRO+feiQXCM7fDgWXZhCl38a4eESkOtAMSXWhTKaWUC1zp6ZcFPgYqAXuA\nLsaYYyJyOzDTGBMmIvWBuWT841IIeN8Y82YO59OevlJK5ZFOrayUUn5Ep1ZWSimVIy36SinlR7To\nK6WUH9Gir5RSfkSLvlJK+REt+kop5Ue06CullB/Roq+UUn5Ei75SSvkRLfpKKeVHtOgrpZQf0aKv\nlFJ+RIu+Ukr5ES36SinlR7ToK6WUH3G66ItIWRFZIyI7RWT1xRWycji2sIgkisgKZ9tTSinlOld6\n+lHAGmNMLWBt5n5OBgLJgK6SkgtZFz32d/pd/E2/i7/pd+E8V4r+Q8C8zO15wMPZHSQidwChwCwg\n16u7+DP9D/pv+l38Tb+Lv+l34TxXiv5txpjDmduHgdtyOO4d4AXgggttKaWUygdFrvamiKwBymfz\n1tCsO8YYIyJXDN2IyIPAr8aYRBEJcSWoUkop1zm9MLqIbAdCjDGHRCQAiDPG1LnsmDHA48A5oBhw\nI7DYGNMrm/PpeL9SSjkhLwuju1L03wB+N8aME5Eo4CZjTI4Xc0WkJfC8MeYfTjWolFLKZa6M6Y8F\n2onITqB15j4icruI/DeHz2hvXimlLHK6p6+UUsr7WH8iV0Q6iMh2EUkVkZds57FFRCqKSJyIbBOR\nrSISaTuTbfpQXwYRuUlEFolIiogki0hT25lsEZGXM/8fSRKRD0XketuZ3EVEYkTksIgkZXkt1w/J\nXmS16ItIYWAy0AEIBLqLSF2bmSxKBwYbY+4EmgL9/Pi7uEgf6sswEfjEGFMXqA+kWM5jhYhUAZ4G\nGhhjgoDCQDebmdxsDhm1Mqu8PCQL2O/pNwF2GWP2GGPSgYVAR8uZrDDGHDLG/JC5nUbG/9i3201l\njz7Ul0FESgPNjTExAMaYc8aYPy3HsuU4GZ2jEiJSBCgBHLAbyX2MMV8Af1z2cq4eks3KdtGvAOzL\nsr8/8zW/ltmjCQY22k1ilT7Ul6Eq8JuIzBGR70VkpoiUsB3KBmPMUWA88DNwEDhmjPncbirrcvuQ\n7CW2i76//2y/goiUAhYBAzN7/H4n60N9+HEvP1MRoAEw1RjTADhBLn7C+yIRqQ4MAqqQ8Su4lIj0\nsBrKg5iMu3KuWVNtF/0DQMUs+xXJ6O37JREpCiwG5htjltnOY9G9wEMishtYALQWkfcsZ7JlP7Df\nGJOQub+IjH8E/FEj4CtjzO/GmHPAEjL+W/Fnh0WkPEDmQ7K/XusDtov+JqCmiFQRkeuArkCs5UxW\niIgAs4FkY8wE23lsMsb8yxhT0RhTlYwLdeuye4rbHxhjDgH7RKRW5kttgW0WI9m0HWgqIsUz/39p\nS8aFfn8WC4RnbocD1+wsXnXunYJmjDknIv2BVWRciZ9tjPHLOxOAZkBPYIuIJGa+9rIx5jOLmTyF\nvw8DDgA+yOwY/Qg8YTmPFcaYzZm/+DaRca3ne2CG3VTuIyILgJZAORHZB7xCxkOxH4tIBLAH6HLN\n8+jDWUop5T9sD+8opZRyIy36SinlR7ToK6WUH9Gir5RSfkSLvlJK+REt+kop5Ue06CullB/Roq+U\nUn7k/wDY8//eOSbGSgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = np.linspace(0, 10)\n", "p = plt.plot(x, f(x), 'k-')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### `quad` 函数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quadrature 积分的原理参见:\n", "\n", "http://en.wikipedia.org/wiki/Numerical_integration#Quadrature_rules_based_on_interpolating_functions\n", "\n", "quad 返回一个 (积分值,误差) 组成的元组:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.integrate import quad\n", "interval = [0, 6.5]\n", "value, max_err = quad(f, *interval)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "积分值:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.28474297234\n" ] } ], "source": [ "print value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "最大误差:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.34181853668e-09\n" ] } ], "source": [ "print max_err" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "积分区间图示,蓝色为正,红色为负:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "integral = 1.284742972\n", "upper bound on error: 2.34e-09\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjXX/x/HXZ86ItFhyk2RpoeSmaCH93KZIIyVSku6y\nVUokJamkkbuivUQphe7QpiS7W52QIpFkxp617Et2M+f6/P6YqSZmzDhn5nzP8nk+Hh6d65xreXce\nfOYz3+u6vpeoKsYYY+JDgusAxhhjwseKvjHGxBEr+sYYE0es6BtjTByxom+MMXHEir4xxsSRkIu+\niCSLyFIRWSEij+TweZKI7BaRhVl/+oR6TGOMMcFJDGVjEfEBrwONgY3A9yIyXlXTjlj1a1VtHsqx\njDHGhC7UTv8yYKWqrlHVdOAD4IYc1pMQj2OMMaYAhFr0KwDrsy1vyHovOwXqi8giEZkkIheEeExj\njDFBCml4h8yCnpcFQEVV3S8iTYFxQLUQj2uMMSYIoRb9jUDFbMsVyez2/6Sqe7K9niwiQ0SktKru\nyL6eiNgkQMYYEwRVzfcQeqjDO/OBqiJSRUROAG4BxmdfQUTKiYhkvb4MkCML/h9U1f6o8uSTTzrP\nECl/7Luw78K+i2P/OV4hdfqqmiEiXYGpgA94R1XTRKRz1udDgZuAe0UkA9gPtAnlmMYYY4IX6vAO\nqjoZmHzEe0OzvR4MDA71OMYYY0Jnd+RGoKSkJNcRIoZ9F3+x7+Iv9l0ET4IZEyoMIqKRksUYY6KF\niKBhPJFrjDEmiljRN8aYOGJF3xhj4kjIV+8YkxtVJS1tG2vXbmHLlq3s2LGTE08sxT/+UZ5//as8\nZcqcQtYtHMaYMLGibwrUqlW/MnjwVCZM+JJVq77E8w4gUo6EhH8gUgrVnXjeb6j+RpEiZfjnP5O5\n4YZk7r77KsqXP9V1fGNinl29Y0KmqrzxxiyefHIQ27bNIDGxCRkZVwFXAeeQ8ySrCqQiMoWEhKkE\nAt9z4YW3M3z4Q9SuXTms+Y2JZsd79Y4VfROSN96YRq9ej7Bv335UuwF3AMF07L/i871KIDCMCy64\nljFj/kOtWlb8jcmLFX0TFl999Qtt2z7I5s0/ofoCmY9RKIjrAnaTmPgqGRmDaNbsST77rAtFitj1\nBsbkxoq+KVSBgMdNN73IuHEDEOmBak+gWCEcaSkJCXdy0kkwffq71K1rs3EbkxO7OcsUmtWrt3LG\nGc0YP34c8AOqfSicgg9wPp43k/37b+Hyy/+P116bUkjHMSa+WNE3+TJs2NdUrVqb7dsvwvP8QJUw\nHDWBQKAbqp/RvXsHWrd+JaipZI0xf7HhHZOndu1G8N57vYGRwDWOUqxBpDnnnluPJUveoEgRn6Mc\nxkQWG9M3BUZVadToP/j975I5g/b5jhPtISGhJVWqlGfZshEkJlrhN8aKvikQ6ekZXHzxvfz88wJU\nJwKnu46UZT8JCddRtWollix5F5/PRihNfAv7iVwRSRaRpSKyQkQeOcZ6l4pIhojcGOoxTeE6fDiD\n6tVvZ8mSNaj6iZyCD1Acz/uCFSvWULv2XXie5zqQMVElpKIvIj7gdSAZuAC4VUSq57LeQGAKOd+e\naSLE4cMBzjuvA7/8sg3PGw+c4jpSDk7C8yawZMlyGjR42HUYY6JKqJ3+ZcBKVV2jqunAB2TepXOk\nbsAnwNYQj2cKUUaGR9WqnVi3biOe9zlwoutIx3Aynvc53377Bbfc8rbrMMZEjVCLfgVgfbblDVnv\n/UlEKpD5g+CNrLds4D4CqSq1a3dhw4Zf8LwvgOKuI+VDaVQn8NFHfXjppa9chzEmKoQ6y2Z+Cvgr\nQG9VVcmcRzfX4Z2UlJQ/XyclJdlzMMOoWbNnWLJkLqpfAye5jnMcqgEf0LNnGy69dBYNGtiduya2\n+f1+/H5/0NuHdPWOiNQDUlQ1OWv5UcBT1YHZ1lnNX4W+DLAfuEtVxx+xL7t6x5GuXd9jyJC+qM4B\nznAdJygJCW9TpMhLbNo0n5Ilo+mHljGhCeslmyKSCCwDGgG/AvOAW1U1LZf1hwNfqOqnOXxmRd+B\nV16ZTo8e/wa+IvNcfPTy+dpTtaqQljbcdRRjwiasl2yqagbQFZgKpAIfqmqaiHQWkc6h7NsUvokT\n03jwwduAj4n2gg8QCLzOsmXf0q3b+66jGBOx7OasOLV27S7OOacugUBvoIPrOAVoEdCYqVPn0KRJ\nVddhjCl0dkeuyVN6eoDy5a9n585z8LxBruMUgiEULTqMbdu+5eSTi7oOY0yhsqmVTZ4aNerLzp37\n8byXXEcpJPeSkVGRa6/9j+sgxkQcK/px5oknxjJ79ig872OgiOs4hUQIBN5k1qyhjB270HUYYyKK\nDe/EkZkzV5GUVA/VScClruOEwUiKFXuZHTvmceKJJ7gOY0yhsOEdk6Pduw9xzTWtgSeIj4IPcAfp\n6WfQosUA10GMiRjW6ceJ6tW7sXz5RjxvLPE1590GoDaTJn1J06Y1XYcxpsBZp2+O8vDDn7Bs2UQ8\n713iq+ADnInIf2jd+h4CAZuG2Rgr+jHu++/X8eKLXVD9ECjpOo4Tqnexf38GHTu+5zqKMc7Z8E4M\nS08PUKZMI/buTcbzeruO49h8RK5j1ao0zjqrlOswxhQYG94xf2re/Hn27lU8zx40Apcg0oLrr+/r\nOogxTlmnH6NGj/6B225rCswHKrmOEyG2Axfw6adTadnyItdhjCkQNg2DYdu2/ZQvX4eMjBSgjes4\nEUVkGKecMpxdu2aT+XgHY6KbDe8YGjbsjerFWME/mmpH9u49SJ8+n7iOYowT1unHmBdf/JKHH74D\n1cWAnbDM2ZckJt7Fzp2pNiGbiXrW6cexjRt/p1evjqi+jRX8Y7kK1ercfvsbea9qTIyxTj+GVKt2\nF6tXC4HAW66jRIFURJJYuXIZZ59tPyBN9Ap7py8iySKyVERWiMgjOXx+g4gsEpGFIvKDiFwV6jHN\n0Z5+ehIrV/6PQOBF11GixAWItKRly6ddBzEmrEJ9Rq6PzGfkNgY2At9zxDNyReQkVd2X9bom8Jmq\nnpvDvqzTD9L69bupUuWfeN5IwH6m5t8moAazZ8/niivOch3GmKCEu9O/DFipqmtUNR34ALgh+wp/\nFPwsJwPbQjymOcJVV/UEmmEF/3idTkLCfdxxRz/XQYwJm1CLfgVgfbblDVnv/Y2ItBCRNGAycH+I\nxzTZPPvsdFatmobnPec6SlTyvIdYvXoi06YtdR3FmLBIDHH7fI3HqOo4YJyINAD+C5yX03opKSl/\nvk5KSiIpKSnEeLHt11/30KfPXagOBU51HSdKlSAh4SE6dUph/foPXIcxJk9+vx+/3x/09qGO6dcD\nUlQ1OWv5UcBT1YHH2GYVcJmqbj/ifRvTP041atzHsmX7CQSGu44S5fYB5zJ27BRuvPFC12GMOS7h\nHtOfD1QVkSoicgJwCzD+iEDnSNb97iJSB+DIgm+O35Ahs0hNHUcgEKsPNw+nkxDpzb332mRsJvaF\nVPRVNQPoCkwFUoEPVTVNRDqLSOes1VoBi0VkIfAqNjdAyHbvPkj37ncCg7CbsAqGame2bl3A++/P\ncx3FmEJlN2dFofr1H2fevKUEAmNdR4kxb1G69Kds3z7FdRBj8s2mYYhxH3+8iG+/fZtA4HXXUWJQ\ne3buTOXdd63bN7HLOv0ocuhQBiVLXs7Bg/cAnVzHiVGDKV16Ktu3j897VWMigHX6MaxVq1c5fPhU\noKPrKDGsEzt2/MCYMQtdBzGmUFinHyVmzvyFhg0vBb4DjprFwhSolylXbjabNtk5ExP57MlZMcjz\nlNNOS2b37qtQPWpOO1Pg9gNnM378/7j++n+6DmPMMdnwTgy6995R/P77ZlQfdB0lThRHpAddu9oM\nnCb2WKcf4ZYt20b16v9EdQJwies4cWQPcDZ+/xwaNqzqOowxubLhnRhTufIdbNhwGp73susocSch\noS/Vqm0mLW2o6yjG5MqGd2LIs89OY/36mXhef9dR4pLndWPp0o9ZvPg311GMKTDW6Ueobdv2c/rp\n/yQQGAw0dR0nbiUkdOPii09i3rwBrqMYkyMb3okRF1/ci0WLNhAIjHYdJc6tAS5mzZrVVK5cwnUY\nY45iwzsxYPTohSxYMJJA4BXXUQxV8PmSufPON10HMaZAWKcfYQ4ezKBkybocOtQNaO86jgHgJ0SS\n2blzNSVKFHMdxpi/sU4/yt1446ukp5cE2rmOYv5Ui4SEi+je/b+ugxgTMuv0I8iXX66mUaPLgLnA\nOa7jmL/5kiJF7uPAgSX4fNYrmchhnX6U8jzlxhs7I/IIVvAj0ZUEAsV45pnJroMYE5KQi76IJIvI\nUhFZIZkV68jPbxORRSLyk4h8IyK1Qj1mLLrzzpHs2bMd1R6uo5gcCZ73EC+8YI+nNNEt1Aej+4Bl\nQGNgI/A9cKuqpmVb53IgVVV3i0gymQ9Sr5fDvuJ2eOfnnzdTq1ZNVKcAdVzHMbk6DJzNhx9OoHXr\ni1yHMQYI//DOZcBKVV2jqunAB8AN2VdQ1W9VdXfW4lzgzBCPGXOaNLkfkQ5YwY90JyDSjYcfftF1\nEGOCFmrRrwCsz7a8Ieu93HQCJoV4zJjyyCPj2LTpRzwvxXUUkw+qd7Nu3UTmz9/gOooxQUkMcft8\nj8eIyJVkPvLpitzWSUlJ+fN1UlISSUlJIUSLfKtX7+T55+9D9QPgRNdxTL6Uwue7nS5dBjFv3kDX\nYUwc8vv9+P1+AoEAEydOPO7tQx3Tr0fmGH1y1vKjgKeqA49YrxbwKZCsqitz2VfcjemfdVZH1q0r\njufZQ86jy2rgMjZvXkvZsie5DmPikKpy1113sWnTJiZOnBjWMf35QFURqSIiJwC3AH97orSIVCKz\n4P87t4Ifj/r3n8batV/iec+6jmKO29n4fA3o2vU910FMnHrhhReYP38+Y8aMOe5tQ745S0SaAq8A\nPuAdVX1WRDoDqOpQERkGtATWZW2SrqqX5bCfuOn0N27cQ6VKNfG8ocA1ruOYoPhJTLyXAweWkJho\nt7uY8Pnss8/o1q0b3377LRUrVrRZNqPBeefdw6pVGQQCw1xHMUFTEhJq06/fAPr0SXYdxsSJH3/8\nkauvvpopU6Zw8cUXA3ZHbsR77rnprFgxiUDALvuLboLndefFF20mVBMeu3fv5qabbmLQoEF/Fvxg\nWKcfRuvX/06VKjXxvLeBJq7jmJAdBCozaZKfpk2ruw5jYpiq0qpVK8qXL8/gwYP/9pkN70Sws8++\ni7VrBc97y3UUU0BE+nLeeVtJS3vDdRQTw15++WVGjx7N7NmzKVq06N8+s6IfoZ54YgpPP90Z1cXA\nqa7jmALzG3ABv/yyiipVSrsOY2LQt99+S4sWLfjuu+8466yzjvrcxvQj0NKl23n66U6ojsAKfqwp\nj8/XjG7d3nUdxMSg33//nbZt2/LWW2/lWPCDYZ1+IfM85YwzWrN1a0U8z2ZojE1z8fnacODASooU\n8bkOY2JIx44dSUxM5K23ch8Stk4/wtxzzyi2bEnF855xHcUUmrqo/oOnnjr+W+KNyc24ceOYOXMm\nL71UsM2idfqF6Ntv11G//sXANKC26zimUL1PiRIj2bVruusgJgZs3ryZCy+8kE8//ZT69esfc13r\n9CPEoUMBmjS5A5EHsYIfD25m9+7FfPFFWt6rGnMMqsqdd95Jp06d8iz4wbCiX0iuvvpZ9u8XVHu5\njmLCoigid9Orl02eZ0IzatQo1q1bx5NPPlko+7fhnULw+utz6NbtRuAHjv14ARNbfgX+ydq1v1Cp\nUgnXYUwU2rp1KzVr1mTChAlccskl+drGrtN3bNWqXVSrVhvPexVo7jqOCTOfrw3Nm1/Op592dx3F\nRKHbbruN8uXL88ILL+R7Gyv6DnmeUr58G7ZtK4vnDXIdxzjxDT5fBw4dWorPZ6OnJv8mTZpE165d\nWbx4MSedlP/nNNiJXIduumkI27Ytw/Oedx3FOFMf1eI895xdxWPyb8+ePdx7770MHTr0uAp+MKzT\nLyAjR86lffvrgTnAua7jGKeGUabMeLZuHZ/3qsYADz74INu3b2fkyJHHvW3YO30RSRaRpSKyQkQe\nyeHz80XkWxE5KCIPhXq8SLRs2XY6dmwNvIUVfANt2bZtDrNm/eI6iIkCixcv5r///S/PPx+eEYJQ\nn5HrA5YBjYGNwPfAraqalm2dfwCVgRbATlXNcSL5aO30Dx/2KFu2GXv2/NOGdcyfEhIeonZtH/Pn\nP+c6iolgqkrDhg1p06YNXbp0CWof4e70LwNWquoaVU0HPgBuyL6Cqm5V1flAeojHikhXXJHCnj37\nbJoF8zeedy8//DCcnTsPuI5iItioUaPYt28fnTt3DtsxQy36FYD12ZY3EEcXpnfr9jE//DASz/sY\nKOI6joko5+LzXUrPnh+4DmIi1O7du+nVqxdDhgzB5wvfRH2hFv3oG48pIGPGLOT117ugOg4o5zqO\niUCBwH2MGvU6nhe3/0zMMTz55JM0a9aMunXrhvW4iSFuvxGomG25IpndflBSUlL+fJ2UlERSUlKw\nuypUS5Zs5t//bgEMwebVMblryuHD3Rg+fC6dOtVzHcZEkNTUVEaNGkVqaupxb+v3+/H7/UEfO9QT\nuYlknshtROY96PM44kRutnVTgD3RfiJ327YDVKrUmIMHr0K1v+s4JuI9T6VKi1m79j3XQUyEUFWa\nNm1KcnIyDzzwQMj7C/sduSLSFHgF8AHvqOqzItIZQFWHisjpZF7VcyrgAXuAC1R17xH7ifiif+hQ\ngDPPvJkdO4rhee9j97aZvG0HziE1dQXVq//DdRgTASZNmkSPHj1YvHgxJ5xwQsj7s2kYConnKdWr\n38/KlUvwvMlA0Ty3MQbA5+vAVVedx7RpvV1HMY6lp6dTs2ZNXnzxRZo1a1Yg+7RpGApJ48bPs2LF\n13jeZ1jBN8cjELiPGTPe5PDhgOsoxrEhQ4ZQpUoVrr32WmcZrOjnw7///Q5+/2BUJwE2Za45XpcA\n5ejXzx6nGM+2b9/O008/zUsvvYRIvhvzAmfDO3m49973efPN3sBXQFXXcUzUeo8SJUaxa9dU10GM\nIw888ADp6ekMHjy4QPdrY/oF6N57P+LNNx8A/gdc4DqOiWoHgUpMnfoNTZpY8xBvVq1aRd26dUlN\nTaVs2bIFum8r+gXknns+ZejQLmQ+1LyW6zgmBoj0pmbNwyxa9JLrKCbMWrduzYUXXsjjjz9e4Pu2\nol8A2rZ9lzFjHgcmYTdfmYKzBriYzZvXUbZs4c6ZbiLH3LlzadWqFcuXL6d48eIFvn+7eidE1133\nAh988BTwNVbwTcGqgs93BQ8+ONp1EBMmqkrPnj156qmnCqXgB8OKfpbDhz3q1OnNpEnvoDoLqOY6\nkolBgUBXPv54sM3HEyc+//xzdu/eTbt27VxH+ZMVfWD9+r2UL38TixbNyir4FfPcxpjgNCY9fT9D\nh85xHcQUsoyMDB599FEGDhwY1lk08xL3RX/GjDWcfXZ9du0qhed9CZRxHcnEtARUu/D00wV72Z6J\nPCNHjqRcuXIkJye7jvI3cX0i9z//mUrfvu1QfRS4H3B3w4SJJ7uAs1i0aCm1atm03LHowIEDVKtW\njY8//ph69Qp3hlU7kZsP27cfpFatHvTteyeqY4DuWME34VOShISbuf/+t1wHMYVk8ODBXHLJJYVe\n8IMRd53+8OGLufvuf+N5VfG8t4DShX5MY472EyLXsnfvLxQvbk9diyW7du2iWrVq+P1+Lrig8G/q\ntE4/F7/9todatR6mY8eryMjonvWIQyv4xpVaiJzNE0987jqIKWDPP/881113XVgKfjBivtNPT1e6\ndv2IYcN6ItKYQGAgULC3QRsTnI84+eQh7Nnjdx3EFJBNmzZRo0YNFi5cSKVKlcJyTOv0sxw+7NG1\n66cUL16bYcOew/PGEAgMxwq+iRwt2bt3BWPHLnYdxBSQZ555hnbt2oWt4AejIJ6clcxfT84apqoD\nc1jnNaApsB9or6oLc1inQDr9LVv28dhjnzBy5It43gl43pPAddiJWhOJRJ7i3HN/ZfnyN11HMSFa\nu3YtderUIS0trcAnVTuWsM69IyI+Mp+R25jMh6R/zxHPyBWRa4GuqnqtiNQFXlXVo05ph1L0Dx3y\nGDLke157bQRr1nyIz1efQKALmT9nrNibSLYJqM4vv/xClSolXYcxIejYsSMVKlSgf//wPjv7eIt+\nYojHuwxYqaprsg7+AXADkP3B6M2BkQCqOldESopIOVXdHOxBVeG77zbzwQezmThxIqtXT0KkJJ53\nG/ATgcCZQf8PGRNep+PzNaV79xF8/nnoD8k2bixdupQvvviCFStWuI6Sp1CLfgVgfbblDUDdfKxz\nJnDMon/w4GFWrdrN8uU7SU1dz7Jla1i9ei1paYvZuXM+qvvw+S4jELgWeBzVc0L8XzHGjUCgGxMn\n3kFGxv0kJsbsabaY1rdvX3r27EnJkpH/21qoRT+/4zFH/uqR43annXYagUCAQ4cOcehQBqolgJJk\nzoVTJetPG+BFTjzxLKePHDOmoHhePcoc/J0pp57CdScU4jX7lSvDokWFt/84tXDhQmbPns3w4cNd\nR8mXUIv+Rv4+O1lFMjv5Y61zZtZ7R+nUqRMigs/no0GDq7nwwitDjGdMNBDG92rO4DHvct2B/YV3\nmGXLCm/fcaxPnz489thjnHRSeJ6R4Pf78fv9QW8f6oncRDJP5DYCfgXmcewTufWAVwr6RK4x0e7g\n7t1ULlWKmaqcV1gHKVoUDh4srL3HpTlz5nDrrbeyfPlyihYt6iRDWK/TV9UMoCswFUgFPlTVNBHp\nLCKds9aZBKwWkZXAUKBLKMc0JhYVK1GCuy+/nNcTbEw/mvTp04e+ffs6K/jBiPk7co2JFhsXLKDm\nxRfzC1CiMA5gnX6B+vLLL+ncuTNpaWkkJoY6Uh48uyPXmChVoU4dmlSowAi7QCHiqSqPP/44/fr1\nc1rwg2FF35gIcn9KCoMAz3UQc0yTJk1i7969tGnTxnWU42ZF35gIcnmnTpQsWpRJroOYXHmeR58+\nfXjqqadIiMJzMNGX2JgYJiJ079CBVyPomarm7z799FN8Ph8tWrRwHSUodiLXmAhzeN8+qpxyClNV\nqVmQO7YTuSELBALUrFmTl156KWKefWsnco2JciecdBJdrrySV6zbjzijR4+mdOnSXHPNNa6jBM06\nfWMi0LalS6lavTrLKMAnQFinH5L09HTOP/983n33XRo2bOg6zp+s0zcmBpQ5/3xuPucc3rTLNyPG\niBEjOPvssyOq4AfDOn1jItSS8eNpfMMNrAEK5H5P6/SDdvDgQapVq8bHH39M3bpHTiTslnX6xsSI\nGs2bU6tUKca4DmIYOnQoF110UcQV/GBY0TcmgvXo1YuXExLyPYe5KXj79u1jwIABYX8iVmGxom9M\nBLumVy8CiYn8z3WQOPbaa6+RlJTEhRde6DpKgbAxfWMi3IjOnRn9zjtMCwRC25GN6R+3Xbt2UbVq\nVb755huqVavmOk6Owvpg9IJkRd+YnB3au5ezTz2VSaqE1Gta0T9uffr04bfffuOdd95xHSVXVvSN\niUEDr72Wn6dN47+hdPtW9I/Lli1bqF69OgsWLKBy5cqu4+TKir4xMWjX2rWcXaUKi/j7s0ePixX9\n49KjRw8yMjIYNGiQ6yjHFLaiLyKlgQ+BysAaoLWq7sphvXeBZsAWVc11KhEr+sYc24MXXUTC4sW8\n4AU58bIV/Xxbu3YtderUITU1lXLlyrmOc0zhvE6/NzBdVasBM7KWczIciIyZiYyJYg+88QbDPY+j\nOitT4Pr160eXLl0ivuAHI5ROfynQUFU3i8jpgF9Vz89l3SrAF9bpGxOa2ytVosaGDfQO5t+Kdfr5\nkpqaSlJSEitWrKBEiUJ5cGWBCmenX05VN2e93gzE3o9EYyJM79df5xVV9rsOEsP69OlDr169oqLg\nB+OYD3cUkenA6Tl89Hj2BVVVEQm5TU9JSfnzdVJSEklJSaHu0piYUqN5c+qVLcu7W7bQ1XWYGDRv\n3jy+//57Ro0a5TpKrvx+P36/P+jtQx3eSVLVTSJSHvjKhneMKXxzR46kdYcOrFDlhOPZ0IZ3jklV\nady4Mbfeeit33nmn6zj5Fs7hnfFAu6zX7YBxIezLGJNPddu1o+qppzLadZAYM3XqVDZu3Ej79u1d\nRylUoRT9AcDVIrIcuCprGRE5Q0Qm/rGSiIwB5gDVRGS9iHQIJbAxBh7r148BIoQ4MYPJEggE6NWr\nFwMGDCAx8Zij3lHPbs4yJgqpKpeffDIP799Pq/xuZMM7uRo5ciRvv/02s2bNQqLswTU2n74xcUBE\neKxXL/onJBDkrVomy4EDB3jiiSd47rnnoq7gB8OKvjFR6vonniChSBE+dx0kyg0aNIhLL72U+vXr\nu44SFja8Y0wUG//YY/QdOJAFnpd3B2fDO0fZvn07559/PrNnz+a8885zHScoNuGaMXFEAwEuKV6c\nPocP0zKvla3oH6V79+5kZGQwePBg11GCZmP6xsQR8flI6dWLfja2f9yWLVvG6NGj/3ZTaDywTt+Y\nKKeex6UnncRjBw9y47FWtE7/b66//noaNmxIz549XUcJiXX6xsQZSUggpU8fUkSs28+n//3vf6Sl\npdGtWzfXUcLOOn1jYoCqUu+UU3hg3z5uzW0l6/SBzBuxateuTUpKCjfeeMzfjaKCdfrGxCERYcBz\nz9FHhMOuw0S4d955h1KlStGyZZ6nvmOSdfrGxJDkMmW4fscO7svp35J1+uzYsYPq1aszefJk6tSp\n4zpOgbBLNo2JYws/+ohrb7mFFcDJR35oRZ+uXbsSCAR44403XEcpMFb0jYlzt555JjV++40+Rz5L\nN86L/o8//sg111xDamoqp512mus4BcbG9I2Jc/1HjuQVz2Ob6yARRFXp2rUr/fv3j6mCHwwr+sbE\nmHMbNaJNjRr08/lcR4kYo0aN4uDBg3Tq1Ml1FOdseMeYGLRt1SouqFqVr1Sp8cebcTq8s2vXLmrU\nqMHYsWM2s6VIAAALl0lEQVSpV6+e6zgFzoZ3jDGUOecc+tx8Mz0SEoj3VurRRx/l+uuvj8mCH4yQ\nOn0RKQ18CFQG1gCtVXXXEetUBN4DygIKvKWqr+WwL+v0jSlA6QcPcmGJEgw4fJjmEJed/uzZs7nl\nlltYsmQJJUuWdB2nUIS70+8NTFfVasCMrOUjpQM9VLUGUA+4T0Sqh3hcY0weihQrxsvPPMNDIhxy\nHcaBQ4cOcffdd/Pqq6/GbMEPRqid/lKgoapuFpHTAb+qnp/HNuOAQao644j3rdM3phBcf/rpNNi6\nlV5FisRVp9+/f3++//57Pv/885h+IlZYr9MXkZ2qWirrtQA7/ljOZf0qwNdADVXde8RnVvSNKQQr\nZs3i8n/9iwUnnEClQ/HR8y9btowrrriChQsXUrFiRddxCtXxFv08H/suItOB03P46PHsC6qqIpJr\n1RaRk4FPgO5HFvw/ZJ/XOikpiaSkpLziGWPyULVBA7onJdFl5ky+UI3prhcgIyOD9u3bk5KSEpMF\n3+/34/f7g96+IIZ3klR1k4iUB77KaXhHRIoAE4DJqvpKLvuyTt+YQnL499+pfd55pLz2GjfffLPr\nOIXq2WefZcaMGUybNo2EhNi/QDHcwzvPAdtVdaCI9AZKqmrvI9YRYGTWej2OsS8r+sYUom+++YbW\nrVvH9JUsixYtonHjxvzwww9UqlTJdZywCHfRLw18BFQi2yWbInIG8LaqNhOR/wNmAj/Bn5cMP6qq\nU47YlxV9YwrZPffcA8Cbb77pOEnBO3ToEJdddhk9evSgffv2ruOEjU24ZozJ1R93p77//vtceeWV\nruMUqEcffZTU1FTGjRsX8+ctsrOib4w5pilTpnD33XezaNEiSpXK9WK7qDJ9+nTat2/PggULKFeu\nnOs4YWVF3xiTp27durF161bGjBkT9V3xr7/+yiWXXMKoUaNi7reX/LC5d4wxeRo4cCA//fQTo0eP\ndh0lJBkZGbRt25Z77rknLgt+MKzTNyZOLVy4kCZNmjB//nwqV67sOk5Q+vbty5w5c5g6dSq+OJ1K\n2oZ3jDH59vzzz/PJJ58wc+ZMihYt6jrOcZkwYQKdO3eOy3H87KzoG2PyTVW5+eabKVmyJG+//XbU\njO//9NNPNGrUiAkTJlC3bl3XcZyyMX1jTL6JCCNGjOC7776Lmmv3N23aRPPmzXn99dfjvuAHwzp9\nYwwrV67kiiuuYOzYsfzf//2f6zi5OnDgAFdeeSVNmzblySefdB0nItjwjjEmKFOmTKFDhw58/fXX\nVKtWzXWco6Snp9O6dWuKFSvG6NGjo2YoqrDZ8I4xJijJycn079+fJk2asH79etdx/iYjI4Pbb7+d\nw4cPM2LECCv4IchzamVjTPy488472b17N1dffTUzZ86kbNmyriPheR4dO3Zk+/btfPHFF1F3lVGk\nsaJvjPmbhx56iF27dpGcnMyMGTOcTtUQCATo3Lkza9euZfLkyRQrVsxZllhhwzvGmKM89dRTXHnl\nlTRo0IB169Y5ybB3715atGjBmjVrmDBhAsWLF3eSI9ZY0TfGHEVEeOGFF+jQoQP169fnxx9/DOvx\nf/31V/71r39RtmxZJk+ezCmnnBLW48cyK/rGmByJCA899BAvvfQSTZo0YfLkyWE57rx587j88stp\n1aoVw4YNo0iRImE5brywSzaNMXmaNWsWbdu2pWXLlgwYMKBQhlrS09Pp378/Q4cOZciQIbRq1arA\njxGLwnbJpoiUFpHpIrJcRKaJyFHPXxORYiIyV0R+FJFUEXk22OMZY9xp0KABP/30E9u3b6d27drM\nnTu3QPe/ePFi6tWrxw8//MCPP/5oBb8QhTK80xuYrqrVgBlZy3+jqgeBK1X1IqAWcGXW4xONMVGm\nVKlSjBo1iqeffpobbriBW2+9lZ9//jmkfaalpdG2bVsaNWrEPffcw4QJEyhfvnwBJTY5CaXoNyfz\ngedk/bdFTiup6v6slycAPmBHCMc0xjh20003sWLFCi666CIaN25Mq1atmD59OocOHcrX9gcOHGDC\nhAm0adOGhg0bUqtWLVatWsVdd91lN12FQdBj+iKyU1VLZb0WYMcfy0eslwAsAM4B3lDVXrnsz8b0\njYky+/btY9iwYXz44YcsWbKEpKQkkpKSqFChAmXLlqVMmTLs3LmT9evXs379eubOncuMGTOoXbs2\nLVu2pGPHjnZlTogKdO4dEZkOnJ7DR48DI7MXeRHZoaqlj7GvEsBUoLeq+nP4XLNPoPTHXx5jTHTY\ntm0b06ZNY86cOWzevJktW7awdetWSpUqRcWKFalYsSI1a9akWbNmnHbaaa7jRi2/34/f7/9zuV+/\nfuGZcE1ElgJJqrpJRMoDX6nq+Xls8wRwQFVfyOEz6/SNMeY4hXPCtfFAu6zX7YBxOYQp88dVPSJy\nInA1sDCEYxpjjAlBKJ1+aeAjoBKwBmitqrtE5AzgbVVtJiK1gBFk/nBJAP6rqs/nsj/r9I0x5jjZ\nfPrGGBNHbD59Y4wxubKib4wxccSKvjHGxBEr+sYYE0es6BtjTByxom+MMXHEir4xxsQRK/rGGBNH\nrOgbY0wcsaJvjDFxxIq+McbEESv6xhgTR6zoG2NMHLGib4wxccSKvjHGxJGgi76IlBaR6SKyXESm\n/fGErFzW9YnIQhH5ItjjGWOMCV0onX5vYLqqVgNmZC3npjuQCthTUvIh+0OP4519F3+x7+Iv9l0E\nL5Si3xwYmfV6JNAip5VE5EzgWmAYkO+nu8Qz+wv9F/su/mLfxV/suwheKEW/nKpuznq9GSiXy3ov\nAw8DXgjHMsYYUwASj/WhiEwHTs/ho8ezL6iqishRQzcich2wRVUXikhSKEGNMcaELugHo4vIUiBJ\nVTeJSHngK1U9/4h1ngFuBzKAYsCpwFhVvSOH/dl4vzHGBOF4HoweStF/DtiuqgNFpDdQUlVzPZkr\nIg2Bnqp6fVAHNMYYE7JQxvQHAFeLyHLgqqxlROQMEZmYyzbWzRtjjENBd/rGGGOij/M7ckUkWUSW\nisgKEXnEdR5XRKSiiHwlIktE5GcRud91Jtfspr5MIlJSRD4RkTQRSRWReq4zuSIij2b9G1ksIqNF\npKjrTOEiIu+KyGYRWZztvXzfJPsHp0VfRHzA60AycAFwq4hUd5nJoXSgh6rWAOoB98Xxd/EHu6kv\n06vAJFWtDtQC0hzncUJEqgB3AXVUtSbgA9q4zBRmw8msldkdz02ygPtO/zJgpaquUdV04APgBseZ\nnFDVTar6Y9brvWT+wz7DbSp37Ka+TCJSAmigqu8CqGqGqu52HMuV38lsjoqLSCJQHNjoNlL4qOos\nYOcRb+frJtnsXBf9CsD6bMsbst6La1kdTW1grtskTtlNfZnOAraKyHARWSAib4tIcdehXFDVHcCL\nwDrgV2CXqv7PbSrn8nuT7J9cF/14/7X9KCJyMvAJ0D2r44872W/qI467/CyJQB1giKrWAfaRj1/h\nY5GInAM8AFQh87fgk0XkNqehIohmXpWTZ011XfQ3AhWzLVcks9uPSyJSBBgLvK+q41zncag+0FxE\nfgHGAFeJyHuOM7myAdigqt9nLX9C5g+BeHQJMEdVt6tqBvApmX9X4tlmETkdIOsm2S15beC66M8H\nqopIFRE5AbgFGO84kxMiIsA7QKqqvuI6j0uq+piqVlTVs8g8UfdlTndxxwNV3QSsF5FqWW81BpY4\njOTSUqCeiJyY9e+lMZkn+uPZeKBd1ut2QJ7N4jHn3ilsqpohIl2BqWSeiX9HVePyygTgCuDfwE8i\nsjDrvUdVdYrDTJEi3ocBuwGjshqjVUAHx3mcUNVFWb/xzSfzXM8C4C23qcJHRMYADYEyIrIe6Evm\nTbEfiUgnYA3QOs/92M1ZxhgTP1wP7xhjjAkjK/rGGBNHrOgbY0wcsaJvjDFxxIq+McbEESv6xhgT\nR6zoG2NMHLGib4wxceT/AZhppPNX9o1QAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print \"integral = {:.9f}\".format(value)\n", "print \"upper bound on error: {:.2e}\".format(max_err)\n", "x = np.linspace(0, 10, 100)\n", "p = plt.plot(x, f(x), 'k-')\n", "x = np.linspace(0, 6.5, 45)\n", "p = plt.fill_between(x, f(x), where=f(x)>0, color=\"blue\")\n", "p = plt.fill_between(x, f(x), where=f(x)<0, color=\"red\", interpolate=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 积分到无穷" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from numpy import inf\n", "interval = [0., inf]\n", "\n", "def g(x):\n", " return np.exp(-x ** 1/2)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "upper bound on error: 7.2e-11\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAADICAYAAAAuo384AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HX504ISYAkhLAIYScgEBBQFsWFXahVwaWi\n8KX2URXrUtf+FLWKfrVKq9alrV8rgnX51gXrF4oKohBRwAqigIAQEWQTFNkNCpl7fn/MEBK2JDDJ\nzSTv5+NxHzP3zpl7P3EehnfOOXOuOecQERERkWPnBV2AiIiISLxToBIRERE5TgpUIiIiIsdJgUpE\nRETkOClQiYiIiBwnBSoRERGR41RioDKzCWa22cyWHKXNE2aWZ2aLzKxrbEsUERERqdxK00M1ERh8\npBfN7GdAG+dcNnAV8FSMahMRERGJCyUGKufcB8C2ozQ5D/hHtO1/gHQzaxib8kREREQqv1jMoWoC\nrCuyvx7IisF5RUREROJCQozOYwftH3I/GzPTPW5EREQkbjjnDs43RxSLHqoNQNMi+1nRY4e48cab\ncM5pi8PtnnvuCbwGbfrsquOmzy9+N3128b2VVSwC1RRgFICZ9QK2O+c2H67h448/xoIFC2JwSRER\nEZHKo8QhPzP7J3AWkGlm64B7gBoAzrmnnXNvmdnPzOxL4AfgV0c+Vx8GDDibLVs2k5AQq9FGERER\nkWCVmGqcc5eWos11pbmY77/Frl31ueiiX/B///ev0rxFKok+ffoEXYIcI3128U2fX/zSZ1e92LGM\nEx7ThcxcZK76TGAA//rX6wwbNqxCri0iIiJSFmaGK8Ok9AACFcDlJCa+wnffbSY1NbVCri8iIiJS\nWmUNVAHdy28CBQV16dt3QDCXFxEREYmhgAKVh++/x8KFn/DYY48FU4KIiIhIjAQ05LffvXjef/PV\nV6to3rx5hdQhIiIiUpI4mUN1gOd1okmT3axdu7pC6hAREREpSZzMoTrA92exfv1Grr/+t0GXIiIi\nInJMAu+hinges8uZN28ePXv2rJB6RERERI4k7ob8Drw+kDp1PuH777/VKuoiIiISqLgb8tvPuTfZ\nvbuAoUMvDLoUERERkTKpNIEKEvH9qbz55r95/fXXgy5GREREpNQqzZDfAVdQo8aLrF+/lgYNGpR7\nXSIiIiIHi9s5VAf4eF42TZo41qz5Es+rRJ1oIiIiUi3E7RyqAzx8/z+sX/8Nw4dfFnQxIiIiIiWq\nhIEKIBPnpvLaa68yfvz4oIsREREROapKOORX1B143h/5/PMltG/fvlzqEhERETlYFZhDVZznnUrt\n2iv47rtNJCYmlkNlIiIiIsVVgTlUxfn+LHbvhjPP7Bt0KSIiIiKHVekDFSTh+x/w8cf/4c477wy6\nGBEREZFDVPohvwPGY3YVM2bMoH///jGrS0RERORgVW4OVXG/IDHx32zYsI7MzMyY1CUiIiJysCoe\nqHw8rzVNm3p89VWeFv0UERGRclHlJqUXF1n0c+3aDVx22cigixEREREB4i5QATTAucm88srLTJw4\nMehiREREROJtyK+o2/C8R1i2bCnt2rWL4XlFRESkuov5kJ+ZDTazL8wsz8xuO8zrmWY2zcw+M7PP\nzezyMtZ8jMYBJ9Oz52ns3bu3Yi4pIiIichhH7aEysxCwAhgAbADmA5c655YXaTMWqOmcG2NmmdH2\nDZ1zBQedK8Y9VAA/4nmN6dSpJZ999kmMzy0iIiLVVax7qHoAXzrn1jjn9gEvA+cf1OYbIDX6PBX4\n/uAwVX6S8P1PWLz4cy666BcVc0kRERGRg5QUqJoA64rsr48eK+oZoKOZbQQWATfErrzSaIlz03n9\n9dcZO/beir20iIiICJBQwuulGaO7A/jMOdfHzFoDM8zsJOfcrkObji3yvE90i4U+wFPce+/VtG9/\nIpdcckmMzisiIiLVQW5uLrm5ucf8/pLmUPUCxjrnBkf3xwC+c25ckTZvAQ845+ZE998DbnPOLTjo\nXOUwh+pgN+N5T/DRR/Po3r17OV9LREREqqpYz6FaAGSbWQszSwQuAaYc1OYLIpPWMbOGQDvgq9KX\nHEuP4txAzjjjLDZu3BhMCSIiIlLtlLgOlZkNAR4DQsCzzrkHzWw0gHPu6eg3+yYCzYgEtAedc/97\nmPNUQA8VRG5P04H09K1s2LCWpKSkCrimiIiIVCVV/F5+pZWP5zWlXbsmfP75Z7rnn4iIiJRJFb+X\nX2ml4Puf8sUXK7nggouCLkZERESquCoaqACa4dy7TJ48mTFj7gi6GBEREanCquiQX1H/AH7FCy88\nz8iRIwO4vkj8mzFjBrNmzWLbtm1ce+215OTksHXrVp577jkyMjJo0qQJAwcODLpMEZGY0Ryqw7oN\ns4eZM+dDTj311IBqEIkfv/3tb3n++ed58MEHOffcc1m0aBHnnHMOq1at4txzz+W///u/mT9/Pvff\nfz8JCQnk5eWxe/duunbtGnTpIiIxUdZAVdLCnlXEOOALzjqrL1999SVZWVlBFyRSaU2ePJk33niD\niRMnMnLkSLKysjj33HMBaN26NZMmTeLkk09mzpw5JCREfoVkZ2czffr0IMsWEQlUFZ5DVZxzbxAO\nt6Fjx85s37496HJEKq1x48bx61//GoD69evTokUL1q07cAeqRYsWccstt3DZZZexadMmAL777rvC\ncCUiUh1VkyG//fbieW2oW/cn1q5dTUpKSsD1iFQuq1atIjs7m48//phTTjml8Pibb77JN998w/bt\n2+ncuTODBg3inXfe4ZFHHqFNmza0bt2am2++OcDKRURiS3OoSpSP57WiUaNEVq/+ksTExKALEqk0\nnnjiCe666y527NiBWal/j4iIVDlah6pEKfj+F2zatJsTT8yhoKAg6IJEKo3c3Fx69uypMCUiUkbV\nMFABpOP7y/j660106XIyvu8HXZBIpfDhhx8WG+oTEZHSqaaBCqARvr+EZcvyOPXU04MuRiRwK1eu\nZMuWLXTp0iXoUkq0ZcsWxo8fz1//+tegSxERAap1oAJojnOfMn/+Qvr3HxR0MSKBmjt3LgAnnXRS\nwJWULDMzk+zsbA3Zi0ilUc0DFUA7nPuIWbNyGTbswqCLEQnM3LlzSUpKom3btkGXElMLFy4MugQR\nqQYUqADognOzmDx5MqNGXR50MSKBmDdvHu3atcPzqtavhXfeeafw+YQJE3j22WcZNmwYixYtCrAq\nEalqtBJfod449zYvvDCYOnXq8Ne/Phl0QSIVZteuXSxbtozhw4cHXUqplWbJl82bN9OoUSMA3n77\nbbp3706nTp3IzMxk1KhRClUiEjNV60/R4zYQeJW//e1vjBlzR9DFiFSYBQsW4Jyjffv2QZdSKlu3\nbmXq1KnMmTOH1atXH7HdlClTGDZsGAB5eXk8/fTTALRp04Y1a9ZURKkiUk2oh+oQFwITeeihy0lN\nTWXMmNuDLkik3H388ccAnHjiiQFXUtzrr7+OmbFw4UJycnJ45513mDBhAhkZGTz88MMlvn/r1q2k\npaUBcM0117B7924A5syZw5AhQ8q19uMxadIkXnzxRRYuXMiWLVto1qwZF1xwAXfccQe1a9cu8f3r\n1q3jpptu4t1338U5x4ABA3jsscdo2rRpBVQvUj2ph+qwRgF/4Y477uDBBx8KuhiRcjd//nwAOnbs\nGHAlByxfvpwzzjiDQYMGMXv2bM4777wyDUmuXr2aVq1aFe4nJCSQnp7O9u3befXVV3nyyZKH9Tdu\n3MiIESPwPI/Fixcf089xLB555BFq1KjBQw89xLRp0/jNb37DU089xcCBA0sc6szPz6dfv36sXLmS\n559/nhdeeIG8vDz69u1Lfn5+Bf0EItWPeqiO6BoA7rjjOnbu3MmDD/4h4HpEys+CBQuoUaNGhX7D\nz/d9nnjiCcLh8CGvtW7dmqFDhwKRuU/9+vUjJSWFQYNKv7zJm2++yZVXXlnsWDgc5v777+eFF16g\nfv36JZ6jcePG3HzzzUybNo3OnTuX+trHa+rUqdSrV69w/8wzzyQjI4Nf/vKX5Obm0rdv3yO+95ln\nnmH16tWsXLmyMFB27tyZ7Oxsnn76aW666aZyr1+kOlKgOqprgFQeeuiX7Nixg7/9TYsIStWzdetW\n1q5dS4cOHQiFQhV2Xc/zuPHGG4/4+pIlS0hKSuLdd9/l7LPPBuDdd99lwIABpTr/nj17qFmzZrFj\n//M//8Ott95Ko0aNeOmllxgxYkSJ55k9ezZnnXVWqa4ZK0XD1H77V7DfuHHjUd87ZcoUTj311GK9\ncy1atKB3795MnjxZgUqknGjIr0QjgX/x1FP/w4gRI4MuRiTmPvvsMwBycnICrqS4t99+m7fffptm\nzZqxZMkSXnzxRXr27Fmq9y5ZsoROnToVO/baa69x++2306lTJ+rXr8+LL75YqnPNnj2bPn36lLX8\nmHv//fcBSvziwNKlSw/7WXbo0IFly5aVS20ioh6qUjofmME//3k227fv4M03/x10QSIxs3/hy/IK\nVCtWrOCFF14gKyuLLVu20KBBA6666qoS3/f//t//O+ZrzpgxgxtuuKHYsYsvvpiLL764xPfOnz+f\niRMn0r59e3zfZ86cOdx3332Fr+/Zs4cnn3ySpKQk5s+fz9VXX81HH33EvHnzuO++++jQocMx130k\nGzZs4O6772bgwIF069btqG23bdtG3bp1DzmekZHBtm3bYl6biEQoUJVaP5yby9tvn84ZZ/Th/fdn\nVrkFEKV62t9DVR5zhD755BNuvvlm3nrrLWrVqsWqVavK/f57zjnC4fAxDV/Onz+fSy+9lLlz59Kg\nQQMmTJiA7/vFeruefPJJrr/+epKTkxk6dChPP/00EyZM4L777mP06NHFAlVBQQHXXHMN+/btK/Ha\nw4cPLxzaLGr37t2cf/75JCYmMnHixDL/TCJSMRSoyqQ7zi1kzpxT6NatBwsXfqxQJXFv8eLFmBld\nu3aN+bl/9atf0bt3b1566SV27dpFw4YN+dOf/hTz6xS1cOFCTj/92G54fsUVV3DVVVfRoEEDINLb\nU3S4zzlH7969SU5OBiK9b48++igJCQns2LHjkPMlJCTw97///ZhqgUhv2LnnnsuaNWt4//33ady4\ncYnvqVu37mF7orZu3UpGRsYx1yIiR2clfQXXzAYDjwEhYLxzbtxh2vQB/gzUALY45/ocpo2Dklc2\njg+r8bxOtGrVlKVLF5GYmBh0QSLHZN++fdSqVYu0tDS+++67mJ575cqVnHjiiWzatKkwoFRmCxYs\noEePHixfvpx27doBcM455zBkyBCuu+66Q9pv2LCBli1bsm3bNmrVqhXzevbt28fQoUP58MMPmTFj\nBj169CjV+/r378/evXv54IMPih3v06cPZsasWbNiXqtIVWRmOOestO2P2kNlZiHgL8AAYAMw38ym\nOOeWF2mTDvwVONs5t97MMo+t9HjSEt9fyVdfdaRVq7asXLmMlJSUoIsSKbMVK1ZQUFDAySefHPNz\nb9++HeCQMLV69WpatmwZ8+sdr1WrVpGWllYYpgoKCvjwww8ZN24cH374YWGvl+/7eJ7He++9x8kn\nn1wYpubMmUPv3r2LnXPfvn1ce+21ZR7y832fESNGkJuby9SpU0sdpgDOO+88br311mL/ndesWcPc\nuXMZN+6Qv4dFJFacc0fcgFOBaUX2bwduP6jNNcB9RztPtJ0DV8W2753nNXCZmY3ctm3bnEi8efXV\nV52ZuTFjxsT83Hv27HH169d3eXl5hccWLFjgxo0bF/NrxcKSJUtcvXr1Cvcfe+wxV6tWLeecc3/4\nwx+cc8699tprrmHDhs4554YNG+ZGjRrlnHNu165d7o9//GPMarn66qudmbm77rrLzZs3r9i2fv36\nwna5ubkuFAq5559/vvDYDz/84Nq0aeM6derkJk+e7CZPnuw6d+7sWrdu7X744YeY1ShS1UUi0tGz\nTdGtpDlUTYB1RfbXAwd/bzkbqGFms4A6wOPOuReOL+bFiwx8fxVbt3agefNWrFixrPBGrCLxYPny\nSGdzecyfSkpK4rXXXuP++++nV69e7N27lyZNmhzXt/fKU05ODjfffDP33XcfqampdO3alb59+/LH\nP/6R0047DYCsrCzOPPNMHnnkEW655RaefPJJnnrqKfLz87n++utjVsu0adMwMx544AEeeOCBYq+N\nHTuWu+++Gyj+B/F+KSkpzJw5k5tuuon/+q//KnbrGfWki5Sfo86hMrMLgcHOuSuj+yOBns6564u0\n+QvQDegPpADzgHOcc3kHncvBPUWO9IluVcFePK8ziYnrmDPngxK/1ixSWVx66aW88sor5OXl0bp1\n66DLEREJTG5uLrm5uYX79957b5nmUJUUqHoBY51zg6P7YwDfFZmYbma3AcnOubHR/fFEhgknHXSu\nKjQp/XB8PG8wMJNXXnmZiy66KOiCRErUpUsX1q1bx/fffx90KSIilUpZJ6WX9J3/BUC2mbUws0Tg\nEmDKQW0mA6ebWcjMUogMCVbD5Xg9fP8dfP9aLr74F8UWAhSpjHzfZ8WKFaVefVxERI7sqHOonHMF\nZnYdMJ3IsgnPOueWm9no6OtPO+e+MLNpwGLAB55xzlXDQLXf40AH7rnnGpYsWcZrr70cdEEih7V6\n9Wp++uknevXqFXQpIiJxr8R1qGJ2oSo/5HewXMzOJienIwsWfKS1qqTSmTx5MsOGDWP69OkMHDgw\n6HJERCqVWA/5yTHrg3NfsHTpGpo0ac6mTZuCLkikmKVLl2JmGvITEYkBBapy1RLfX8u2bXVo0aJV\n4U1oRSqDRYsWkZOTQ2pqatCliIjEPQWqclebcPgL9u49g+7de/Dqq68GXZAIAJ9++ilnnXVW0GWI\niFQJClQVwsO56fj+9VxyyXDGjr036IKkmtu+fTurVq3izDPPDLoUEZEqQYGqQv0Z+Dv33nsfw4Zd\niO/7QRck1dTs2bMxM/r37x90KSIiVYK+5ReI2ZidTVZWIxYs+M8hN48VKW+jR49m5cqVzJo1q/DY\nQw89RNu2bVm4cCGjRo2ibdu2AVYoIhIsfcsvLpyJc+vYsMEjK6sZ06dPD7ogqeJ83+fss8/mgw8+\nYPfu3UyaNInRo0cXvj5nzhxWrlzJBRdcwG9+8xt+97vfBVitiEj8UaAKTCa+n8e+fRcxePAQbr1V\n/4BJ+dm9eze5ubls3LiR22+/ndatWzN8+PDC12fNmkWPHj0AaNKkCfPnzw+qVBGRuKRAFSgPeBF4\njkcf/TNdu55Cfn5+0EVJFZSamsoDDzzAVVddxcKFC5k0qditNtm8eTMpKSmF+6FQiO3bt1d0mSIi\ncUuBqlIYhXPLWbz4axo2PIHPPvss6IKkCrr11lvZsWMHc+fOpVmzZsVe832fUChUuF9QUFBsX0RE\njk6BqtLIxve/IT//ZLp1O5nHH3886IKkGmnSpAk//PBD4X44HKZOnToBViQiEl8UqCqVBHx/Js7d\nx4033syQIedoaQWpEAMGDCjsGc3Ly6N79+4BVyQiEl+0bEKlNQfPO5vMzFTmz//okCEakVi74447\n6NSpE59++ilXXnkl2dnZQZckIhKYsi6boEBVqe3E807DbCXPPTeBkSNHBl2QiIhItaBAVSX9FvgL\nffv2Z9q0N0lMTAy6IBERkSpNgarK+g+eN4SkpALefHMKffr0CbogERGRKksrpVdZPfH9b9mzpx99\n+/bj8st/pQnrIiIilYR6qOLSG5hdRv36dcnNfY/27dsHXZCIiEiVoh6qamEYzm1my5amdOyYwz33\njA26IBERkWpNPVRx73HMbqFt23bMnj2LBg0aBF2QiIhI3FMPVbVzA859RV7ejzRunMWzzz4bdEEi\nIiLVjnqoqpRbgD/Tq9dpvP32VNLT04MuSEREJC6ph6paewRYwPz5a8jMbMCjjz4adEEiIiLVgnqo\nqqw7MRtHixatmD79Td1GREREpAzUQyVRD+DcWr7+OoV27U5k9OjfaN0qERGRclJioDKzwWb2hZnl\nmdltR2nX3cwKzOyC2JYox64xvv8Zzk1g/PjnyMioz8yZM4MuSkREpMo5aqAysxDwF2Aw0AG41MwO\nWUUy2m4cMA0odfeYVJRf4vvb2LnzVPr3H8CgQUPIz88PuigREZEqo6Qeqh7Al865Nc65fcDLwPmH\naXc9MAn4Lsb1Scwk4dxUIJf33ltA3br1tMSCiIhIjJQUqJoA64rsr48eK2RmTYiErKeihzTzvFI7\nE9/fzN69V3DFFVfRoUNnVq1aFXRRIiIicS2hhNdLE44eA253zjkzM4465De2yPM+0U0qngc8CdzI\nihXnk52dzdChF/Dii8+TkpISdHEiIiIVLjc3l9zc3GN+/1GXTTCzXsBY59zg6P4YwHfOjSvS5isO\nhKhMIB+40jk35aBzadmESusNPO8KPG83v//9ndx9991BFyQiIhKosi6bUFKgSgBWAP2BjcDHwKXO\nueVHaD8R+Ldz7l+HeU2BqlLzgfsw+wNpaak899yznH/+4abLiYiIVH0xXYfKOVcAXAdMB5YBrzjn\nlpvZaDMbfXylSuXiAWNxbivbt5/J0KHDaN++EytWrAi6MBERkUpPK6XLEeTheRfi+59z7rnn8b//\n+yK1a9cOuigREZEKoZXSJUay8f3FwP/x1ltzSU/P4M4779Rq6yIiIoehQCUlOI9w+FvC4bt48MGH\nSU2ty+OPPx50USIiIpWKhvykDH4ErsfsOVJTU3n44XFcccUVQRclIiIScxryk3KUBDyDczvYseMc\nrrrqaurVa8hLL70UdGEiIiKBUqCSY5ACPI9zW9i69QxGjhxFw4ZNeOONN4IuTEREJBAKVHIc0onc\nwnEz333XlQsuuJCsrBZMnz496MJEREQqlAKVxEBm9MbL69m4MZvBg4fQsmU2s2fPDrowERGRCqFA\nJTHUGOdmAKv5+usTOOusPjRt2oJXX3016MJERETKlQKVlIPmODcbWMOGDR245JJLqVs3k0cffVTr\nWImISJWkZROkAuwkstzCP6lZM5HrrvsNDzzwAImJiUEXJiIiclgxvTlyLClQCewFfo/n/RWzvQwf\nPpy//OUJ0tPTgy5MRESkGK1DJZVYIjAO399JOPwnXn55GhkZ9Rg0aDBr164NujgREZFjpkAlAfCA\nGwiHv8W5V5g5cwXNm7cgJ+ck3nrrraCLExERKTMFKgnYRYTDq4E5LFtWm3PO+TlpafUYM2YMP/74\nY9DFiYiIlIrmUEklsxO4Hc97AdhD3779ePzxP9OxY8egCxMRkWpEc6gkzqUCf8P3d+H7/yA3dw05\nOZ1o1qwV48eP17ILIiJSKamHSuLAKuAGzKaTmFiDyy4bzsMPP0xGRkbQhYmISBWlHiqpgloDU3Fu\nDz/9dBvPPz+VevUyadeug3qtRESkUlAPlcSp/2D2eyCXUMgYMKA/Dz74B7p06RJ0YSIiUgWoh0qq\niZ449w7O/UhBwZ+ZMeNLunbtRr16Dfnd737H7t27gy5QRESqEfVQSRWyCfg9odBrhMM7ycnpzF13\njeGSSy4JujAREYkz6qGSaqwR8Azh8HZgJkuXpnHppSOoWTOZwYOHMHv27KALFBGRKko9VFLFFQBP\n4Hnj8f0vSEpKYdCgAfz+93dxyimnBF2ciIhUUro5ssgR5QN/JhT6B+Hwl9SqVYdzzhnC3Xf/XguH\niohIMeUy5Gdmg83sCzPLM7PbDvP6CDNbZGaLzWyOmXUuS9EiFSMFuJNweCWwnR9++C2vv/4ROTk5\npKXVY9SoX7Jq1aqgixQRkThUYg+VmYWAFcAAYAMwH7jUObe8SJtTgWXOuR1mNhgY65zrddB51EMl\nldQW4EFCoVcIhzeQmprB2Wf355ZbbqFnz55BFyciIgGI+ZBfNCzd45wbHN2/HcA599AR2tcFljjn\nsg46rkAlcWAT8Aih0L8Ih1dTs2YSvXr14pprruaiiy7C8/Q9DhGR6qA8hvyaAOuK7K+PHjuSXwNv\nlbYAkcqlEfAnwuFVQD4//fQHPvhgO8OHj6BGjZp06tSFRx55hPz8/KALFRGRSqQ0PVQXAoOdc1dG\n90cCPZ1z1x+mbV/gr0Bv59y2g15zcE+RI32im0g88IE3MHsKs3n4/h6ysppz0UVDue6662jdunXQ\nBYqIyHHIzc0lNze3cP/ee++N+ZBfLyJzovYP+Y0BfOfcuIPadQb+RSR8fXmY82jIT6qQBUSGBmcS\nDn9LUlItunc/hZEjL2PUqFEkJSUFXaCIiByH8phDlUBkUnp/YCPwMYdOSm8GzARGOuc+OsJ5FKik\nitoJPIPZK5gtwfd/olGjJgwa1J/rrruW7t27B12giIiUUbmsQ2VmQ4DHgBDwrHPuQTMbDeCce9rM\nxgPDgLXRt+xzzvU46BwKVFJNLAT+Rig0g3B4HQkJieTk5PCLX1zIlVdeSWZmZtAFiohICbSwp0il\nshd4CbMX8LwFhMO7SEmpQ5cuJzF06Hn8+te/JiMjI+giRUTkIApUIpXat8AEzP6N5y0mHN5NrVqp\ndOlyEhdcMJTLL79cAUtEpBJQoBKJK5uAZzGbiuctIRz+gVq10ujWrQvnn38uI0aMoFGjRkEXKSJS\n7ShQicS1jRzowVpOOLyLxMRk2rRpQ79+ZzF8+HBOPfVULTAqIlLOFKhEqpSdwMvAv0lIWEhBwSbM\noEGDE+jVqzvnn38eF198MbVr1w66UBGRKkWBSqRK84HZwCt43gfAKnz/R1JSUsnObkPv3r0YOnQo\nffv2JSEhIeBaRUTilwKVSLWzHvgn8B4JCZ8TDm/COZ86ddI58cS2nHXWGQwdOlRDhSIiZaBAJSLA\nUmAS8D6h0FLC4S2YOdLSMsnJac9pp/XinHPO4fTTT1fIEhE5DAUqETkMH/iESMj6kISEPMLh73HO\np1atVJo3b0737t3o168fP//5z7V0g4hUewpUIlIGK4ApwGxCoWU4txHf/5GEhJo0bNiIzp07cNpp\npzFw4EC6d++u3iwRqTYUqETkOO0E3gbexfM+wWwN4fAOwCcpqRYnnNCYnJz29OzZg0GDBnHyyScr\naIlIlaPuEvbiAAAKXElEQVRAJSLlZBWRoDUPz/scs3WFQSs5uXZh0OrWrSunn346vXv3JikpKeCa\nRUSOjQKViFSwPGA6MBfP+xzP20g4vB3nwiQkJJKeXo8WLZqSk9OBHj160K9fP7Kzs9WrJSKVmgKV\niFQSW4GZwDxgMaHQV8BmwuF8AJKTa5GZ2YCWLZvSvv2JdO3ald69e9OhQweFLREJnAKViFRyPrAc\nmAUsBPJISFiH72/B9/MBR40aSaSn16Vp0ya0a5dNp06d6NmzJz169NCq8CJSIRSoRCTOrQU+ILLM\nwzJCoTXAZnx/F86FMfNITq5F3boZNG3ahNatW9KxY0e6dOlCz549teSDiMSEApWIVGE/EglaC4El\nwCpCoXXAlmjgKsDMo2bNZNLS0mnYsAEtWjSlVatWtG/fnpNOOomTTjpJk+VFpEQKVCJSjRUAn7G/\ndwu+wvPW43nf4dx2fH8PzvnR0JVCWloajRo1pGnTxrRo0YI2bdpw4oknctJJJ9GoUaNgfxQRCZQC\nlYjIUeUDnwKLiczlWoXnfYPnbcG5Hfh+Ps4VABAK1SA5uRZpaWk0bFifrKzGNG/evDB8tW/fntat\nW2sSvUgVpEAlInLcfGAdkWHFSOiCr/G8TdHgtQvn9uD7+wCHmUeNGjVJSalNenoaDRpk0rhxI7Ky\nsmjWrBktW7akTZs2tG3blpSUlCB/MBEpJQUqEZEKtZ1I6FpJJHitBTbieZvxvG0HhS8fMEKhBBIT\nk0hJqUV6ehqZmRk0atSARo0accIJJ5CVlUWLFi1o1aoVTZs2JSEhIcCfT6R6UqASEam0CoCviYSv\n1dHnG4DNeN4WPG8b8EM0gP1UOPR4IITVJDk5hTp16lC3bjoZGenUq5dB/fr1OeGEEwrDWLNmzWje\nvLl6w0SOgwKViEiVkk8kfK0m0vu1HvgG+BbYRii0A7NdQD7O/Yjv740GscjvW88LEQrVIDGxJklJ\nydSqlUKdOrWpWzeN9PQ06tWrR7169WjYsCENGzakcePGNG7cmCZNmpCamqr5YVJtKVCJiAiR4cWN\nwBoi88E2R7ctwPfAdjwvEsbMIr1izv2Ec/uKBTKwaChLoEaNRGrWPBDMUlNrk5pah7S0VNLS0khP\nT6du3brUq1ePzMxMMjMzadiwIY0aNSIzM1PhTOKKApWIiMTITiK9Yd8QCWPfAt8RCWTbiMwf24nn\n7cYsH7M9wE9Fgtk+nPOJhLsIM68woCUkJJCQsD+kJZGcnETt2inUrl2L2rVrkZqaSu3atalTpw5p\naWmFW0ZGBnXr1iUjI4OMjAwyMzNJTEys6P84UsUpUImISCW0l0gg20QkkG0lEsq2AjuIhLNdRELc\nD5j9gOflY/YT8BOwF+f2AvtwLoxzBdGwVvTfFcPM8LxIaPO8/aGtBomJNUhMjIS35OSkwi0lJZnk\n5GRSUlIKt9q1axcGuTp16hQGu/2P+3vjUlJS1OtWhcU8UJnZYOAxIASMd86NO0ybJ4AhRAb7L3fO\nfXqYNgpUcS0X6BNwDXJsctFnF89y0edXEp9IIPs++ridSEjbv+0Cdke3H6KP+cCPmO3ffsJsL5Hg\ntxcoiIa2AiCMc2HAJ/Jv5qH/lpl5mFlhD5zneTjnSEysWaQ3LoEaNYqGu0QSE2uQlFQz+vxAb93+\nx0jPXXLhY9HnSUlJpKSkkJycTK1atQoDYdHnCnzHrqyB6qjfxTWzEPAXYACRr6LMN7MpzrnlRdr8\nDGjjnMs2s57AU0CvY6peKrFc9Es9XuWizy6e5aLPryQekBHdysa5yHZsCoj0qO2ILo+xE9hJOLw/\nuE1i794+RMJbJMDBniKPe6PP9wL5mG3D8/YB+zArYH+wgzD7A14kPIYLh1IP9NIdOewdEOnB2/94\nIAB6hY8HthChkEcoFBmeDYU8EhISCIVC0XAYokaNGtHH/WGxeHAsuoVCocLnCQkJJCYmRsNlYuF+\n0ef72xZ9XrNmzWLPi55j/2PRLSEhoUIDZUmLm/QAvnTOrQEws5eB84ksurLfecA/AJxz/zGzdDNr\n6JzbXA71ioiIVBIJHD3IfQ38rtRncw7C4RiUVcgnEtj2B7o9OLenyGPRgPcj+4dWiz8vuu076LEg\n+rygyPYjZgWYhaOhMBIGzfzo8zD7Q2HksXhAjGwHB8TIdvj9/cqSiiOdTmZFn+8PmnAgeJZNSYGq\nCZGvh+y3HuhZijZZRGYwFpOaem6ZC5TK4ccfV5CU9EnQZcgx0GcX3/T5xa/q8dmFolvNoAs5RCSk\nFUSHaw/08EWGciPhLjIfL1xsH3zC4e+JhMbSKylQlTbyHRzlDvu+nTunlvJ0Uhnt3ZsXdAlyjPTZ\nxTd9fvFLn131UVKg2gA0LbLflEgP1NHaZEWPFVOWiV0iIiIi8aSk2VoLgGwza2FmicAlwJSD2kwB\nRgGYWS9gu+ZPiYiISHVy1B4q51yBmV0HTCcySPqsc265mY2Ovv60c+4tM/uZmX1J5GsNvyr3qkVE\nREQqkQpb2FNERESkqir3BRrMbLCZfWFmeWZ2W3lfT2LHzJqa2SwzW2pmn5vZb4OuScrOzEJm9qmZ\n/TvoWqT0okvQTDKz5Wa2LDqlQuKEmY2J/u5cYmb/a2aV72twUsjMJpjZZjNbUuRYhpnNMLOVZvaO\nmaUf7RzlGqiKLAw6GOgAXGpm7cvzmhJT+4CbnHMdiSzWeq0+v7h0A7AM3aog3jwOvOWcaw90pvj6\nf1KJmVkL4Eqgm3OuE5EpM8ODrElKNJFIVinqdmCGc64t8F50/4jKu4eqcGFQ59w+YP/CoBIHnHOb\nnHOfRZ/vJvILvXGwVUlZmFkW8DNgPIcubyKVlJmlAWc45yZAZD6rc25HwGVJ6e0k8gdpipklACkc\n5tvvUnk45z4gcnPJogoXLo8+Dj3aOco7UB1u0c8m5XxNKQfRv7i6Av8JthIpoz8TWarZD7oQKZOW\nwHdmNtHMFprZM2aWEnRRUjrOua3AI8BaYCORb7+/G2xVcgyK3vVlM9DwaI3LO1BpiKEKMLPawCTg\nhmhPlcQBM/s58G30ZuXqnYovCUA34G/OuW5EvkF91OEGqTzMrDVwI9CCSK9+bTMbEWhRclxcyTdK\nLPdAVZqFQaUSM7MawOvAi865/wu6HimT04DzzGw18E+gn5k9H3BNUjrrgfXOufnR/UlEApbEh1OA\nuc65713kfif/IvL/o8SXzWbWCMDMTgC+PVrj8g5UpVkYVCopi9wd8llgmXPusaDrkbJxzt3hnGvq\nnGtJZELsTOfcqKDrkpI55zYB68ysbfTQAGBpgCVJ2XwB9DKz5Ojv0QFEvhgi8WUK8Mvo818CR+1U\nKOnWM8flSAuDluc1JaZ6AyOBxWb2afTYGOfctABrkmOnIfj4cj3wUvSP0VVo0eS44ZxbFO0NXkBk\n/uJC4O/BViVHY2b/BM4CMs1sHXA38BDwqpn9GlgD/OKo59DCniIiIiLHp9wX9hQRERGp6hSoRERE\nRI6TApWIiIjIcVKgEhERETlOClQiIiIix0mBSkREROQ4KVCJiIiIHKf/D3hzxzX6sL0TAAAAAElF\nTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "value, max_err = quad(g, *interval)\n", "x = np.linspace(0, 10, 50)\n", "fig = plt.figure(figsize=(10,3))\n", "p = plt.plot(x, g(x), 'k-')\n", "p = plt.fill_between(x, g(x))\n", "plt.annotate(r\"$\\int_0^{\\infty}e^{-x^1/2}dx = $\" + \"{}\".format(value), (4, 0.6),\n", " fontsize=16)\n", "print \"upper bound on error: {:.1e}\".format(max_err)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 双重积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "假设我们要进行如下的积分:\n", "\n", "$$ I_n = \\int \\limits_0^{\\infty} \\int \\limits_1^{\\infty} \\frac{e^{-xt}}{t^n}dt dx = \\frac{1}{n}$$" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def h(x, t, n):\n", " \"\"\"core function, takes x, t, n\"\"\"\n", " return np.exp(-x * t) / (t ** n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "一种方式是调用两次 `quad` 函数,不过这里 `quad` 的返回值不能向量化,所以使用了修饰符 `vectorize` 将其向量化:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from numpy import vectorize\n", "@vectorize\n", "def int_h_dx(t, n):\n", " \"\"\"Time integrand of h(x).\"\"\"\n", " return quad(h, 0, np.inf, args=(t, n))[0]" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [], "source": [ "@vectorize\n", "def I_n(n):\n", " return quad(int_h_dx, 1, np.inf, args=(n))" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(array([ 1.97, 1. , 0.5 , 0.2 ]),\n", " array([ 9.804e-13, 1.110e-14, 5.551e-15, 2.220e-15]))" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "I_n([0.5, 1.0, 2.0, 5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "或者直接调用 `dblquad` 函数,并将积分参数传入,传入方式有多种,后传入的先进行积分:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.integrate import dblquad\n", "@vectorize\n", "def I(n):\n", " \"\"\"Same as I_n, but using the built-in dblquad\"\"\"\n", " x_lower = 0\n", " x_upper = np.inf\n", " return dblquad(h,\n", " lambda t_lower: 1, lambda t_upper: np.inf,\n", " x_lower, x_upper, args=(n,))" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(array([ 1.97, 1. , 0.5 , 0.2 ]),\n", " array([ 9.804e-13, 1.110e-14, 5.551e-15, 2.220e-15]))" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "I_n([0.5, 1.0, 2.0, 5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 采样点积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### trapz 方法 和 simps 方法" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy.integrate import trapz, simps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`sin` 函数, `100` 个采样点和 `5` 个采样点:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_s = np.linspace(0, np.pi, 5)\n", "y_s = np.sin(x_s)\n", "x = np.linspace(0, np.pi, 100)\n", "y = np.sin(x)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVfXa//H3VwbJEXFONCewHHFi2iioleijZc5alkpP\n46mec9kvw9JMq2OWdTyp2DH1yk6KRT5qZWIqVGqS5aygR0RFHFLxoCDz/v7+EHmMQKa9WXu4X9fl\ndfaGtdf+tDzeLO691vdWWmuEEEI4llpGBxBCCGF5UtyFEMIBSXEXQggHJMVdCCEckBR3IYRwQFLc\nhRDCAZVb3JVSK5RSF5VSh+6wzT+UUv9WSh1QSvW0bEQhhBCVVZEz95VAeFnfVEoNBTpqrX2Ap4Ao\nC2UTQghRReUWd631T8DVO2zyEPBp0bYJgKdSqrll4gkhhKgKS/TcWwGptz0/C3hbYL9CCCGqyFIf\nqKoSz2VNAyGEMJCrBfaRBrS+7bl30df+QCklBV8IIapAa13yBLpcljhz3wg8DqCUCgT+o7W+WNqG\nWmu7/fPGG28YnsFR81+6dIlZs2b94XlcXFyp2dPT0+nRowcdO3bk1VdfJTQ0lOHDh1OvXj1iYmJK\n3f/Jkyc5fPhw8fPly5cTHx8vx17y28Wfqir3zF0ptQYIBZoopVKBNwC3omL9sdZ6k1JqqFLqBJAF\nTKlyGuE0zpw5Q+vWrVFK0bBhQ1q2bInWGqUUTZo0ISws7E+vOXHiBAMHDqRp06ZMmDABFxcXAHr3\n7k3Dhg154oknSElJ4eWXX/7D69q1a/eH576+vjRr1qz4eWJiIp06daJWLbntQziOilwtM0FrfbfW\n2l1r3VprvaKoqH982zZ/0Vp31Fr30FrvtW5k4QgmT57MiRMnAHBzc+OZZ55BqbJ/8/zpp5/w9/fH\nx8eH4cOHFxf2tm3bAtCxY0cmTZrE3LlzefbZZzGbzWXuKyQkBF9fX+Dmb5MvvvgiaWl/6iQKYdfk\nVKWCSjuTtCdG59+3bx8//vhj8fNt27bh4+NTodcWFBQwZMgQwsLC6N+//x9+CNx+Vt6iRQsiIiJY\nv349Q4cOJS8vr9x9K6X4/vvvad365sdGFy5cYOnSpRX9z6oQo499dUl++6Sq09Op1BsppWvqvYTt\n2b59OxkZGTzyyCOVet1bb73FvHnzGD169J/aK2XJyckhJiaGOnXqEBcXh5eXV4XfLyUlhbi4OKZO\nnVqpnEJYi1IKXYUPVKW4C6vIy8vjww8/ZNq0abi6Vv6iLLPZTEREBBs2bGDChAl/6JFXRGFhIZs2\nbeLixYts27atuA1TWQsXLmTIkCFVfr0Q1VXV4i5tGWEVbm5uuLi4kJ2dXenX5uTkcP/997Nlyxam\nTp1a6cIO4OLiwrBhw+jUqRMBAQH88MMPld4HQNOmTSt15i+ErZAzd2ExFy5c4NixY4SGhlZ5Hxcv\nXmTAgAForRk5ciTu7u7VznXw4EFiY2OJiopi0qRJVd7PyZMnSUpKYujQodXOJERFyZm7MFxqaiq7\nd++u8uuPHDlCz549adCgAWPHjrVIYQfo3r07Y8aM4dlnn2XOnDlV3s+VK1e4eLHUWziEsDly5i6q\n5caNG7i4uFC7du1q7WfLli2MGTOGoKAgAgMD73hZZFVdunSJ1atXM3z4cFauXFnt69rPnj2Lt7cs\noySsS87chSFmzpzJF198Ua19LFu2jJEjRxIeHk5QUJBVCjvc7J9HRESwbds2Bg4cWKXPA265cOEC\no0aNorCw0IIJhbAcOXMX1ZKbm4u7u3uVC3JkZCSLFi1i3LhxxdeaW1teXh7r1q0DIC4ujhYtWlRp\nP2azWe5qFVYnZ+6ixsyfP5/k5GQAateuXaXCXlhYyLhx41i2bBlTpkypscIO4O7uztixY/H09KRn\nz54cPny4Svu5Vdjz8vJ48sknuXr1TmMPhKhZUtxFpfn4+FCvXr0qvz4zM5N+/fqRkJDA1KlTady4\nsQXTVUytWrUYPHgwfn5+mEwmYmNjq7wvNzc3hg8fjqenpwUTClE90pYRFZKdnc1dd91V7f2kpqYy\ncOBAPDw8ePjhh3Fzc7NAuuo5evQo3377LR988AFPPfVUtfeXnp4u18YLi5G2jLCqhx9+mH379lVr\nH3v37qV37940b96ckSNH2kRhB+jcuTPjx49n2rRpTJ8+vVr7KigoYNCgQXLJpDCcnLmLCsnIyKBh\nw4ZVfv3XX3/No48+Sv/+/enbt68Fk1lOeno6q1evJiwsjOjo6OKVJysrLy/PYtfoCyFn7sLi9u/f\nX3y5YHUK+0cffcSECRMYPny4zRZ2AC8vL6ZOncqePXsICQkhMzOzSvu5Vdi11nz99dd3XH5YCGuR\n4i7KtGLFCg4cOFCtfbz00ku8/vrrPProo3ax+FadOnV47LHHyMrKws/Pj9TU1PJfVIbc3FzWr19f\n5R8SQlSHtGWEVeTn5zNq1CgSEhKYOHGi3V1JYjab2b59O0lJSWzevJnevXsbHUk4KWnLCIv44Ycf\nqn22npGRQWBgIEeOHGHKlCl2V9jh5qWS999/P4GBgQwYMIANGzZUa39Xr17l7bfflhaNqDFS3MUf\nXLp0qVo345w8eZIePXpgNpuZOHGiRS6fNFKfPn0YPnw4EydO5O9//3uV9+Pm5kbjxo2ttrSCECVJ\nW0ZYzK5duxg+fDjdu3cnNDTUoQrZuXPniI6O5vHHH2fhwoWy7ICoMdKWEVWWlJTEkiVLqrWPtWvX\nMnjwYPr3709YWJhDFXaAu+++m6lTp/Lll1/y0EMPkZ+fX+V97dixgwULFlgwnRB/JsVdULduXZo3\nb17l18+bN4+IiAhGjhyJn5+fBZPZFk9PT6ZMmUJSUhIBAQFVbl917NiRwMBAC6cT4o+kLSOqzGw2\n89RTTxETE8PEiROr9QPCnhQUFPDtt99y5coVtm/fTocOHYyOJByYtGVEpeTm5jJ9+nSysrKq9Pqc\nnBwGDx7Mpk2biIiIcJrCDuDq6spDDz1Ehw4d6Nu3Lzt37qzyvmbPnk1CQoIF0wlxkxR3J+Xq6kqn\nTp2qdDXL5cuX6du3L2fOnGHy5Mk0aNDACgltm1KK0NBQQkNDGTx4MNHR0VXaz4MPPmgXN3cJ+yNt\nGVEpiYmJPPDAA9x9992Eh4dXef0VR3Ly5Em++uorIiMjee2114yOIxyMtGVEhaxdu7bKQ6zj4uII\nCgqic+fODB06VAp7kfbt2/PEE0/w/vvvM3Xq1CrdqHTlyhWeffbZal2FI8TtpLg7mQYNGlC/fv1K\nv27lypUMHz6cBx54AJPJ5HCXOlZXs2bNiIiIIDY2lgceeICcnJxKvb5Ro0Y88MADuLq6WimhcDbS\nlhHlmjVrFh9++CFjxozhnnvuMTqOTcvNzWXdunW4uLgQHx9P06ZNjY4k7Jy0ZUSZsrOz+fjjj6ns\nD9dbSwgsXryYyZMnS2GvgNq1azN27Fjq1atHjx49SExMrPQ+Nm7cyLZt26yQTjgTKe5OIDMzk8uX\nL1fqNVlZWYSGhrJjxw4iIiJo0qSJldI5HhcXF4YMGUK3bt0ICgpi69atlXp9o0aN7HKxNWFbpC0j\n/uTcuXMMGDAANzc3RowYYTPj8OzRkSNH2LRpE//4xz+YOnWq0XGEHZK2jPiTvXv38vvvv1fqNfv3\n76dnz540adKEUaNGSWGvpi5dujBu3LjioSWVkZuby6efflrpdpoQIMXdoW3fvr1SQ62//fZb+vfv\nT58+fXjwwQdl5UMLadOmDVOmTGHp0qWMHz+ewsLCCr2uoKCAQ4cOkZeXZ+WEwhFJW0YAsGTJEl55\n5RWGDx/Ovffea3Qch5SVlcXatWvx9vZmy5Yt1K1b1+hIwg5IW0YUO3v2bKW2nzZtGq+++ioTJ06U\nwm5FdevWZdKkSWRkZNCzZ0/S0tIq/Nrk5GSuXbtmxXTC0ZRb3JVS4UqpJKXUv5VS00v5fhOl1Gal\n1H6l1GGl1GSrJBUVcu3aNR566CGys7PL3bawsJARI0bw2WefMXXqVFq1alUDCZ2bm5sbo0aNolmz\nZvTq1avCbbN//vOf/Pzzz1ZOJxzJHdsySikX4BhwP5AG7AEmaK0Tb9tmNlBbax2plGpStH1zrXVB\niX1JW6aGmM3mcvvl165dY9CgQVy5coWxY8fa/Tg8e7Rnzx5+/PFHVq9ezbBhw4yOI2yUtdoy/sAJ\nrfUprXU+EA08XGKb88CtZQEbAFdKFnZhfVrr4jVNyivsp06dws/Pj/z8fB599FEp7Abp27cvw4YN\nY/z48Xz00UcVft2NGzesmEo4ivKKeysg9bbnZ4u+drtlQBel1DngAPCS5eKJivriiy94+eWXy90u\nISGBPn360KZNG0aMGCFrmRisU6dOTJw4kddee43/+Z//KXf7H3/8kUcffbQGkgl7V96/7Ir0UWYA\n+7XWYUqpDsD3SqkeWuvrJTecPXt28eOwsDDCwsIqEVXcyahRoxgwYMAdt/nqq6+YPHkyAwcOpFev\nXjWUTJSnVatWTJ06ldWrV5OSkkJMTEyZ9xf069eP3r1713BCUZPi4+OJj4+v9n7K67kHArO11uFF\nzyMBs9b63du22QS8rbXeWfR8GzBda/1riX1Jz91A77//PrNnz+aRRx6hY8eORscRpcjOzubLL7/E\ny8uLrVu30rBhQ6MjCRtgrZ77r4CPUqqtUsodGAdsLLFNEjc/cEUp1RzoBJysbBBRNfPnz2fv3r1l\nft9sNvPMM8/w1ltvMWnSJCnsNuyuu+5i4sSJFBQU4OfnR0pKSpnb5uTkMGLECLk8UpTpjsW96IPR\nvwCxwFFgrdY6USn1tFLq6aLN3gH6KKUOAFuBV7TW6dYMLf5Pz549adOmTanfy8vLY+jQoWzYsIGI\niAhatGhRw+lEZbm6uvLwww9zzz330KdPnzIHq3h4ePD//t//q9La/MI5yB2qDurKlSsMHDiQ7Oxs\nRo0ahYeHh9GRRCXt3buXbdu2sXLlSsaMGWN0HGEQuUPVyWzatKnMNUqOHz+On58fbm5ujB8/Xgq7\nnerVqxcjR45kypQpzJ8/v8ztYmJiOH/+fA0mE/ZAirsdys3NJTo6utTrnX/44QcCAgLo1KkTw4YN\nkzmndq5Dhw48/vjjvPPOOzz11FOlzme9dOkSV65cMSCdsGXSlnEgn332Gc8++yyDBw+me/fuRscR\nFnTt2jXWrFlD165d+eabb3B3dzc6kqghVW3LSHG3M1lZWaWuJjhnzhzmz5/PmDFjaNu2bc0HE1aX\nm5tLTEwMtWvXJi4u7k/TsXJycrhx4wZeXl4GJRTWID13J3Ds2DEefPDBPwxvMJvNPP744yxcuJDJ\nkydLYXdgtWvXLv4Mxc/Pj6SkpD98f+nSpSxfvtygdMLWyJm7nbn9zD07O5vw8HBOnjzJ+PHjqVev\nnsHpRE3QWrNz50727NnD+vXri+/0LiwspFatWihV6ZM8YcOkLeNkLly4wIABA1BK8cgjj0gP1gkd\nOnSIzZs3s2jRIiZPnmx0HGEl0pZxYDt37mTBggXFzw8ePEjPnj3x9PRkzJgxUtidVLdu3Rg7diwv\nvPDCH9Zt2rVrFy+9JOv3OTtZEtAOtG/fvvhxbGwsY8eOJTg4mMDAQANTCVtwzz33MHnyZD766CNO\nnDjBqlWr6NGjBw0aNCj/xcKhSVvGjixdupSXX36ZYcOGcd999xkdR9iQzMxM1q5dS9u2bYmNjZU1\n+h2ItGUcUEpKSvGczenTp/PKK68wYcIEKeziT+rVq8ekSZO4cuUKfn5+nDt3juvXr7NxY8l1/oSz\nkOJuw3766SfWr1/P6NGjWbFiBVOmTMHb29voWMJGubu7M3r0aBo3bkyvXr349ddf2b59O/Ibs3OS\ntowNy8zM5P777+fixYuMHTuWOnXqGB1J2ImEhAR27txJdHQ0Q4YMMTqOqAZpyziQuLg4zpw5Q48e\nPbhx4waPPvqoFHZRKQEBAQwZMoQxY8awdOlScnJyLDLdR9gPOXO3Qd26deP06dP4+fkxYMCAcgde\nC1GWs2fPEh0djbu7O8899xxz5841OpKoJDlzdxDbtm3j2LFjmEwmBg0aJIVdVIu3tzcREREopVi/\nfr3RcUQNkjN3G3FrKO7ixYu5fPkyoaGhALRt25Z27doZnE7Yq5SUFE6dOkV+fj67du0iIiICb29v\nGVBvR6p65i43MdmIsLAwtm/fTm5uLv369WPAgAFGRxIOoF27dsUnB2fPnuX7778nOTkZV1f5p+/o\n5Hd+G/Lpp5/So0cPacUIq/D29ub8+fPs37/f6CiiBkgVsRE7duwgPT2d/v37y7K9wip8fHzw9/fn\nzTffNDqKqAHSc7cBN27cYMiQIbi5udGvXz+j4wgHdv36dZYsWUJcXBwBAQFGxxEVIFfL2LERI0aQ\nkJBA3759jY4iHFz9+vXp0qUL4eHhXL161eg4woqkuNsApRR9+/bFw8PD6CjCCQQHB5Obm0teXp7R\nUYQVSXE3WHJyMj/99JMs3ytqjJeXF76+vsycOdPoKMKKpLgb6OOPP+b555+nR48eMiJP1CiTycSa\nNWt47LHHuH79utFxhBVIcTdQfn4+P/zwg5y1ixrXrFkzWrduTWZmpkzyclBS3A20d+9eOnfujKen\np9FRhBMymUzEx8djNpuNjiKsQIq7Aa5cuUJ6ejpffPEFwcHBRscRTsrb25smTZowf/58jh07ZnQc\nYWFynbsBHnvsMTIzM0lOTmb06NFGxxFO7OTJk3z33Xd07NiRTZs2yW+RNkjWlrEjixcvpk2bNowf\nP97oKMLJtWvXDg8PD0aPHi2F3cFIW8YA77//Ps2bN+fuu+82OopwckopTCYTH374ofTeHYwU9xq0\nfPlydu3aRVRUFCEhIUbHEQIAX19f8vPz+fjjj5kwYQI5OTlGRxIWIMW9BrVs2ZKNGzdSr1497rnn\nHqPjCAFArVq1MJlMfPDBB0yePBk3NzejIwkLkOJeg8LDw1m1apWctQub06VLF9LT07l+/TouLi5G\nxxEWIMW9BqSnp1NYWMiyZcsA6Nixo8GJhPgjFxcXgoODmTt3LlprEhMTjY4kqqnc4q6UCldKJSml\n/q2Uml7GNmFKqX1KqcNKqXiLp7Rz7733HitWrODdd9/FZDKhVKWvahLC6vz8/EhJSeF///d/efrp\np8nPzzc6kqiGO17nrpRyAY4B9wNpwB5ggtY68bZtPIGdwGCt9VmlVBOt9eVS9uW017lrrVmzZg0v\nvvgizz//vExaEjZr165dZGRksHv3bqOjiCLWWs/dHzihtT6ltc4HooGHS2wzEfhKa30WoLTC7uyU\nUrzzzjuYTCYp7MKm9e7dm0OHDpGQkGB0FFFN5VWaVkDqbc/PFn3tdj6Al1IqTin1q1JqkiUD2rMt\nW7bw7bff8t1333H27Fm6d+9udCQh7qh27dr4+/szY8YMTp8+zQsvvGB0JFFF5RX3ivRR3IBewFBg\nMDBTKeVT3WCOwNPTE09PT2bNmkVQUJBMnBd2wd/fn59//pmrV68SHh6Os7ZT7V151SYNaH3b89bc\nPHu/XSpwWWudDWQrpX4EegD/Lrmz2bNnFz8OCwsjLCys8ontiL+/Pzt27CAxMZGXXnrJ6DhCVEid\nOnXo1asXs2bNYuPGjUbHcTrx8fHEx8dXez/lfaDqys0PVAcB54Bf+PMHqvcCi7h51l4bSADGaa2P\nltiX03ygmpubC9z8FTc0NBRXV1f69+9vcCohKu7WIO2jR4/Spk0bzp49S5s2bYyO5ZSs8oGq1roA\n+AsQCxwF1mqtE5VSTyulni7aJgnYDBzkZmFfVrKwO5uNGzfy4osvcvDgQfbs2SODr4XdqV+/Pl27\ndmXGjBns2LGDGTNmGB1JVJIs+Wslubm5jBw5kuvXrzNo0CCj4whRaenp6XzyySckJyfTvHlzuT/D\nINa6FFJU0dmzZ4mLiyMgIMDoKEJUiZeXFz4+PrzxxhtS2O2QFHcLunTpElFRUQDMmDGD7t27y+Br\nYdeCg4NZvXo1GRkZxMfH8+mnnxodSVSQFHcLysrKwtXVlXPnzvH1118TFBRkdCQhqqV58+a0bt2a\nuXPn0qJFCzp06GB0JFFB0nO3goiICPbt28fDD5e8mVcI+3P27FliYmI4f/48Hh4eRsdxOtJzN1hB\nQQEAV69eZe3atTL4WjgMb29vGjduzPz58wHIy8sjKyvL4FSiPFLcLeDy5cv4+flRWFjIm2++Sdu2\nbWnatKnRsYSwGJPJxKJFi8jPz+f1119n7dq1RkcS5ZC2jIWkp6dz11130bJlS8aNGyfzUYVD0Vqz\nfPlypk2bxvPPPy/TmmqQtGUM5uXlxd/+9jeaNWsmhV04HKUUISEhLFiwQCY12Qkp7tUUExNDVlYW\n+fn5REVFYTKZjI4khFX4+vqSl5dXfDnksmXLSEtLMziVKIsU92owm83ExcWhtWbhwoXUrVtXBl8L\nh1WrVi1CQkL429/+BoCbmxs5OTkGpxJlkZ67BZjNZlq3bk1YWBi+vr5GxxHCagoLC1m0aBHLli1j\n1KhRRsdxCtJzN9Ann3yC1hofH1nGXjg2FxcXTCYTc+bMKf5aYWGhgYlEWaS4V9GTTz7Jzp07MZvN\nMvhaOJUePXqQkpLC1q1bKSgooHv37qSnpxsdS5QgbZkqOnXqFM2aNeObb77hueeek8HXwqns3LmT\n69ev8/PPP3P58mWaNGlidCSHVdW2jBT3aurWrRvt27enV69eRkcRosbk5uaycOFCtm/fLiufWpn0\n3GvI6dOn+f333wHYvHkzqampMvhaOJ1bg7Rfe+01ADIzM9m0aZPBqcTtpLhX0rZt24iJiQGQwdfC\nqfn7+7Nr1y6OHDlCTk4OGzdulGHaNkTaMlW0c+dOBg8ezEsvvYS7u7vRcYQwxJYtW2jSpAkbNmww\nOorDkp57DQsLC8PFxUUGXwundu3aNaKiojh69Cht27Y1Oo5Dkp67lR0+fJjZs2cDcPDgQX755RcZ\nfC2cXoMGDejatWtx7z02Nrb4sTCWnLlX0O+//86RI0cYMGAAw4YNIyMjg/vvv9/oWEIY7tYg7VOn\nTqGU4tq1a7Rv397oWA5D2jI15OTJk3Tt2pXnn39e5qMKUWTdunUEBQWxdOlSo6M4HGnLWNHVq1eL\nH0dGRsrgayFKMJlMfP7551y7dg24OSxe7lo1lhT3cly7do2AgAByc3O5cOGCDL4WohTNmzfH29ub\nuXPnAvDee+/xww8/GJzKuUlbpgLy8/Nxc3PjySefZO/evTL4WohSpKamsm7dOs6dOyeDtC1I2jJW\n5Obmxn/+8x+io6Nl8LUQZWjdujWNGjXivffeMzqKQIr7Hf3rX/8qnjQjg6+FKF9ISAiLFi2isLAQ\nrTWRkZFkZGQYHcspSXG/g3PnzuHq6kp2djYrV66Us3YhytGuXTvc3d1ZtGgRSim6du2K2Ww2OpZT\nkp57BbzxxhusXr2axx57zOgoQti8xMREdu3axalTp2QZbAuQnrsF3f5DKD8/nyVLlhASEmJgIiHs\nR6dOncjNzWXVqlXFX8vLyzMwkXOS4l6KefPmFd+MsXDhQurUqSODr4WooFq1amEymYoHaaempuLv\n7y8rRtYwacuUIjMzk+zsbBo3bkybNm0IDQ2VwddCVMKtQdqffPIJI0eO5Nq1azRo0MDoWHZJ2jIW\nVK9ePZo2bcry5cspLCyUwddCVJKLiwvBwcHFg7SlsNc8Ke63ycnJYf/+/cDNvvu8efNk8LUQVeTn\n58fJkyfZtm0bACkpKezevdvgVM5Divttjh8/zuLFiwH48ssvycjIoEuXLganEsI+ubm5ERQUxMyZ\nMwFITk7m0KFDBqdyHtJzL0O3bt1o164dvXv3NjqKEHbr1iDtuLg4/P39jY5jl6zWc1dKhSulkpRS\n/1ZKTb/Ddn2VUgVKqZGVDWFrYmNjOXPmDD169DA6ihB2rXbt2vTt21cGeBjgjsVdKeUCLALCgc7A\nBKXUfWVs9y6wGbC7BrXWmhdeeKF4aV8ZfC2E5QQEBLBz504SExMBeOWVV9i5c6fBqRxfeWfu/sAJ\nrfUprXU+EA2UtiTiC0AMcMnC+WqE2WwmODgYT09Pfv75Z44cOSLtGCEspE6dOvTs2ZPIyEgAJk6c\nSLdu3QxO5fjKK+6tgNTbnp8t+loxpVQrbhb8qKIv2U9jvYiLiwsTJkxAKcWMGTPw9/fH3d3d6FhC\nOIzAwEC2bNnC6dOn8fPzk0sja0B5xb0ihfrvwKtFn5Yq7Kwtc+PGjeLHhw4dIiEhQT74EcLCGjRo\nQJcuXf7Qe7+14qqwjvKaymlA69uet+bm2fvtegPRRdeCNwGGKKXytdYbS+5s9uzZxY/DwsIICwur\nfGILe/HFFwkPD2f06NFERkbSu3dv7rrrLqNjCeFwgoOD+eSTT7h06RJ169YlPDycX375Rf69lRAf\nH098fHy193PHSyGVUq7AMWAQcA74BZigtU4sY/uVwNda63WlfM8mL4UsKCjAbDaTlpZGly5deO65\n56hfv77RsYRwSOvWrSM4OJioqCi01nKDYAVY5VJIrXUB8BcgFjgKrNVaJyqlnlZKPV21qLbF1dUV\nd3d3IiMj6datmxR2IawoODi4eJC2FHbrctqbmNLS0khKSmLQoEFcuHCBDh068N///d80atTI6GhC\nOLTo6GiGDx/O/PnzSUpK4sCBA4wbN87oWDZLFg6rpHPnzhXfCj1z5kx8fX2lsAtRA0wmE5988gk5\nOTm4uLiQm5trdCSH5LRn7rdkZGTQqlUrnnjiCZo1a2Z0HCGcwqpVq5gyZQqvv/660VFsnpy5V9Gb\nb77JPffcI4VdiBpkMpmKB2nfYosnf/bM6Yp7Tk4O//Vf/0V2djbZ2dmsWLECk8lkdCwhnEr79u1x\nc3NjyZIlAEybNo2YmBiDUzkWp2vLmM1mdu/eTXBwsAy+FsJAiYmJ/Pzzz6SkpHD+/HmaNWuGm5ub\n0bFsjrRlKqhWrVoEBweTn59PVFSUnLULYZBOnTqRk5PDZ599RqtWraSwW5hTFfeLFy9iNpsB+Oij\nj7jrrrsLC4QnAAAQlUlEQVRo27atsaGEcFK3Bmm/8847wM2e+2+//WZwKsfhVMV91qxZfPPNN5jN\nZj744AMZoSeEwbp27cqlS5dYv349+fn5vPrqq1y/ft3oWA7BqXrut95/+fLlzJw5k6efflqKuxAG\n27NnD2lpaezbt8/oKDZJeu4VoJRCKSWDr4WwIX5+fiQnJxMXF2d0FIfiFMU9OTmZL7/8Erg5+Prq\n1at07tzZ4FRCCLg5SDswMLD4hqb9+/cTFRVVzqtEeZyiuN+4cYOcnBwA5syZg8lkwsXFxeBUQohb\n+vTpw/79+/n1119p3Lgx3t7eRkeye07Vc9+yZQtjxozhxRdflPmoQtiYuLg43N3d2bJli9FRbIr0\n3Ctg5syZMvhaCBvl7+/Pjh07SEpKAm7ecHj78gSichy6uOfm5hIYGEhmZia7d+/myJEj9OrVy+hY\nQohS1K1bl549e/Lqq68C8NRTT7Fx458GuokKcvi2zLFjx+jUqRMDBw4EIDQ0tMYzCCEq5tq1a0RF\nRXHs2DHq1auHp6en01/VJm2ZMnTq1InDhw+ze/du+vbta3QcIcQd3BqkPWPGDBo1auT0hb06HLa4\nnz59muzsbAAiIyPp1asXderUMTiVEKI8QUFBrFu3jsuXL1NQUMD27duNjmSXHLa4L1++nI0bN5KS\nksK2bdsIDAw0OpIQogIaN25Mx44dmTVrFgUFBSxZskSmNVWBw/fcJ06cSHJyMkOHDq3x9xZCVM2F\nCxf4/PPPOXfuHPXq1TM6jqGk516KixcvsmHDBoKCgoyOIoSohBYtWnD33Xczd+5co6PYLYcr7qmp\nqSxatAiQwddC2LOQkBCWLVtGbm4u8fHxxVObRMU4XHE3m814eXmRkZHBmjVrCA4ONjqSEKIKWrdu\njaenJwsWLKBt27Zyj0olOWzPfdq0aXz33XeMGzeuxt5TCGFZycnJbNmyhbS0NKddD0p67vzfeu3Z\n2dksX75cRugJYefat2+Pq6tr8SqRubm5xdPUxJ05THHXWhMUFERaWhrvvvsuTZo0kZXlhLBzSilC\nQkJ4//33MZvNjBgxgl9++cXoWHbBodoyqamptGjRglatWjF06FDatWtn1fcTQlif2WwmKiqKBQsW\nMHLkSOrWrWt0pBolbRlufgCzaNEiPDw8ZPC1EA7i1iDtt99+2+kKe3U4RHE/c+YMly9fLh58HRIS\nImtSCOFAunXrVnzfSlZWFrGxsUZHsnkOsbB5bGwsubm51KlTh/z8fHx8fIyOJISwIBcXF4KDg3nz\nzTcJCQnhiy++4MEHH5STuDtwqJ67j48P3bt3p3v37lZ9HyFEzcvPz2fhwoV8/fXXhIWFGR2nxjh9\nzz0mJob09HS6dOlidBQhhBW4ubkRFBRUPEhb3JldF/esrCxee+01tNbMmTOH4OBgp73RQQhn0Lt3\nb/bt28dvv/3G559/ztq1a42OZLPsurjn5eXh4+PD1q1bOXXqFH5+fkZHEkJYkYeHB3379iUyMpIe\nPXrIkgR34BA998DAQBo2bCjryAjhBLKysli0aBF79+7l3nvvNTqO1Tldz/3WLcgJCQkcPnyY3r17\nG5xICFET6tati5+fHzNmzABuFnvxZxUq7kqpcKVUklLq30qp6aV8/1Gl1AGl1EGl1E6llNUvVxkz\nZgy7du1ixowZ9O3bl9q1a1v7LYUQNiIwMJDNmzdz+vRpevfuzcWLF42OZHPKbcsopVyAY8D9QBqw\nB5igtU68bZsg4KjWOkMpFQ7M1loHltiPRdsyV65c4cyZM5hMJl544QWZjyqEk/n666/p3Lkz//zn\nP/Hw8DA6jtVYsy3jD5zQWp/SWucD0cDDt2+gtf5Za51R9DQBsPqKXY0bN+aNN96QwddCOKng4GC+\n+uorMjMzjY5ikypS3FsBqbc9P1v0tbJEAJuqE+pO0tPTOXHiBCkpKWzdulUGXwvhpBo3bkyHDh14\n4403uHDhAj/++KPRkWxKRZYfqHAvRSk1AJgKlLqQ+uzZs4sfh4WFVekuswMHDrBp0ybS0tLo2rUr\n9evXr/Q+hBCOwWQy8dlnnzF+/Hh27NhB//79jY5UbfHx8cTHx1d7PxXpuQdys4ceXvQ8EjBrrd8t\nsV13YB0QrrU+Ucp+LNZzv3jxIu3bt+fJJ5/Ey8vLIvsUQtinNWvWMGLECObNm2d0FKuwZs/9V8BH\nKdVWKeUOjAM2lnjzNtws7I+VVtgtbdasWfj4+EhhF0JgMplYtmwZeXl5RkexKeUWd611AfAXIBY4\nCqzVWicqpZ5WSj1dtNksoBEQpZTap5Sy+KgUrTWvv/46p0+fZvXq1TJCTwgBQJs2bWjYsCHvv/8+\ns2fP5tdffzU6kk2wmztUCwoKWLZsGSdOnJDB10KIPzhx4gRbt27lyy+/pFOnTjRt2tToSBbj8Heo\nurq6MnnyZFasWCFn7UKIP+jQoQMuLi4cOHDAoQp7ddhFcb+11MB7771H48aNZfC1EOIPlFKYTCbe\ne+89zGazLEmAnRT3N998k8WLF7No0SI5axdClOree+/lxo0bLF68mL59+xafFDoruyjukZGR/Oc/\n/8HDw4N27doZHUcIYYNuDdJesmQJe/fupVYtuyhvVmMX//Xu7u4sXboUk8kkMxOFEGXq1q0bFy5c\n4Pvvvzc6iuFsekC21prffvuNgwcPkpeXh6+vr9GRhBA27NYg7dmzZ9O+fXsApx29adNn7mlpabz1\n1lvMmzePkJAQp/81SwhRvp49e3L8+HGio6NJSkoyOo5hbLpaent7M2nSJC5fvuy0P32FEJXj5uZG\nYGAgcXFxjBo1yug4hrHp4g4wZ84cTCaTDL4WQlRYnz592LdvH/v27TM6imFstuf+8ccf4+LiQkpK\nCsOHDzc6jhDCjnh4eNCnTx+mT59O8+bNiYqKol69ekbHqlE2e+bu5eVFVFQUQUFBuLm5GR1HCGFn\nAgIC2LFjB6Ghobi62ux5rNXYbHFv06YNSUlJMvhaCFEltwZpb9q0yaHH8JXF5hYOu7XNAw88QGFh\nYZUGegghBEBGRgZLly7l+PHjeHl52eVITodZOGzLli088sgj7Nq1C39/f6PjCCHsWMOGDencuTOP\nP/44zzzzjNFxapTNnbmbzWaGDBlCZmYmDz74YA0kE0I4sitXrrB8+XKSk5Np0aKF0XEqzWHO3M+c\nOcNPP/0kg6+FEBZxa5D2W2+9ZXSUGmVTxf3AgQNERkbStWtXGjRoYHQcIYSDuDVI+5tvviEjI8Po\nODXCZoq71pq//vWvbNiwgaCgIKPjCCEcSIsWLWjZsiXvvPMOqampRsepETZT3JVS+Pr64uvrK4Ov\nhRAWZzKZOHbsmNMsQGgzxf3atWt8/vnnBAcHGx1FCOGA2rRpQ4MGDfjggw+MjlIjbKK4x8XFMXny\nZLy9vWnevLnRcYQQDiokJIQPPviAv/71r0ZHsTqbKO6urq58//33MkJPCGFVHTp0wNXVlRs3bhgd\nxepsorjHxcXRrFkzWrdubXQUIYQDU0rRr18/YmNjHX7GquHFPT8/n8WLFxMSEmJ0FCGEE7g1SHv1\n6tXk5eUZHcdqDC3uGRkZtG3bFjc3Nxl8LYSoEbcGab/00kv861//MjqO1Rha3OvXr4/Wmn79+sng\nayFEjenWrRuFhYU0bdrU6ChWY2hxX7VqFYWFhU5z3akQwjbcPkjbURlW3FNSUnj77bcxmUwy+FoI\nUeN69erF8ePHWbBggdFRrMKwqjp9+nTOnz9P165djYoghHBibm5u+Pv7M2/ePHJycoyOY3GGFffj\nx48TFhYmg6+FEIbx9/fnxo0bJCYmGh3F4gwp7lu3buXkyZP4+fkZ8fZCCAH83yDtyMhIo6NYXI0X\n98uXL/Pkk08SGBgog6+FEIYLCAggLi6O+fPnGx3Fomq8uCckJHD+/Hn69OlT028thBB/UrduXe67\n7z42b95sdBSLqvHi/ve//52goCBq165d028thBClGjhwILt37yYtLc3oKBZTo8U9MTGRnTt3EhAQ\nUJNvK4QQd9SwYUPuu+8+h+q9l1vclVLhSqkkpdS/lVLTy9jmH0XfP6CU6lnWvgYNGoSvry916tSp\nTmYhhLC44OBgPv/8c/bu3Wt0FIu4Y3FXSrkAi4BwoDMwQSl1X4lthgIdtdY+wFNAVFn7u3r1KmFh\nYdXNbIiUlBSjI1SLPee35+wg+Y1W0fxNmjTh3nvvZeXKlVZOVDPKO3P3B05orU9prfOBaODhEts8\nBHwKoLVOADyVUqVO3OjatSuNGjWqZmRjnDp1yugI1WLP+e05O0h+o1Umf79+/Vi1ahXfffed9QLV\nkPKKeyvg9mmyZ4u+Vt423qXtTEboCSFsWcuWLWnWrJlD9N5dy/m+ruB+Si7pWOrrWrZsWcHd2R5X\nV1e7vsLHnvPbc3aQ/EarbP7OnTvz7bffUlBQgKtreSXSdimty67fSqlAYLbWOrzoeSRg1lq/e9s2\nS4F4rXV00fMkIFRrfbHEvir6g0IIIcRttNaVXhO9vB9LvwI+Sqm2wDlgHDChxDYbgb8A0UU/DP5T\nsrBXNZwQQoiquWNx11oXKKX+AsQCLsByrXWiUurpou9/rLXepJQaqpQ6AWQBU6yeWgghxB3dsS0j\nhBDCPln8DlVL3vRkhPLyK6XClFIZSql9RX9eNyJnaZRSK5RSF5VSh+6wjU0e+/Ky2/JxB1BKtVZK\nxSmljiilDiulXixjO1s9/uXmt+W/A6WUh1IqQSm1Xyl1VCn1tzK2s9XjX27+Sh9/rbXF/nCzdXMC\naAu4AfuB+0psMxTYVPQ4ANhtyQw1kD8M2Gh01jLy9wN6AofK+L4tH/vystvscS/K1wLwK3pcDzhm\nZ//fr0h+W/87qFP0v67AbiDEXo5/BfNX6vhb+szdojc9GaAi+eHPl37aBK31T8DVO2xis8e+AtnB\nRo87gNb6gtZ6f9HjTCARuLvEZrZ8/CuSH2z77+BG0UN3bp6opZfYxGaPP1QoP1Ti+Fu6uFv0picD\nVCS/BoKLfq3bpJTqXGPpqs+Wj3157Oa4F11d1hNIKPEtuzj+d8hv038HSqlaSqn9wEUgTmt9tMQm\nNn38K5C/Usff0lfoW/SmJwNUJMdeoLXW+oZSagiwHvC1biyLstVjXx67OO5KqXpADPBS0RnwnzYp\n8dymjn85+W3670BrbQb8lFINgVilVJjWOr7EZjZ7/CuQv1LH39Jn7mlA69uet+bmT8c7beNd9DVb\nUG5+rfX1W78+aa2/A9yUUl41F7FabPnY35E9HHellBvwFfAvrfX6Ujax6eNfXn57+DsA0FpnAN8C\nJScC2fTxv6Ws/JU9/pYu7sU3PSml3Ll509PGEttsBB6H4jtgS73pySDl5ldKNVdKqaLH/ty8nLS0\n3pgtsuVjf0e2ftyLsi0Hjmqt/17GZjZ7/CuS35b/DpRSTZRSnkWP7wIeAPaV2MyWj3+5+St7/C3a\nltF2ftNTRfIDo4FnlVIFwA1gvGGBS1BKrQFCgSZKqVTgDW5e9WPzx7687NjwcS9iAh4DDiqlbv2j\nnAG0Ads//lQgP7b9d9AS+FQpVYubJ62faa232UvtoQL5qeTxl5uYhBDCAdX4DFUhhBDWJ8VdCCEc\nkBR3IYRwQFLchRDCAUlxF0IIByTFXQghHJAUdyGEcEBS3IUQwgH9f9NqCj+GhaK3AAAAAElFTkSu\nQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = plt.plot(x, y, 'k:')\n", "p = plt.plot(x_s, y_s, 'k+-')\n", "p = plt.fill_between(x_s, y_s, color=\"gray\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "采用 [trapezoidal 方法](https://en.wikipedia.org/wiki/Trapezoidal_rule) 和 [simpson 方法](https://en.wikipedia.org/wiki/Simpson%27s_rule) 对这些采样点进行积分(函数积分为 2):" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Trapezoidal Integration over 5 points : 1.896\n", "Simpson Integration over 5 points : 2.005\n", "Trapezoidal Integration over 100 points : 2.000\n" ] } ], "source": [ "result_s = trapz(y_s, x_s)\n", "result_s_s = simps(y_s, x_s)\n", "result = trapz(y, x)\n", "print \"Trapezoidal Integration over 5 points : {:.3f}\".format(result_s)\n", "print \"Simpson Integration over 5 points : {:.3f}\".format(result_s_s)\n", "print \"Trapezoidal Integration over 100 points : {:.3f}\".format(result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 使用 ufunc 进行积分" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Numpy` 中有很多 `ufunc` 对象:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "numpy.ufunc" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(np.add)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "accumulate(array, axis=0, dtype=None, out=None)\n", "\n", "Accumulate the result of applying the operator to all elements.\n", "\n", "For a one-dimensional array, accumulate produces results equivalent to::\n", "\n", " r = np.empty(len(A))\n", " t = op.identity # op = the ufunc being applied to A's elements\n", " for i in range(len(A)):\n", " t = op(t, A[i])\n", " r[i] = t\n", " return r\n", "\n", "For example, add.accumulate() is equivalent to np.cumsum().\n", "\n", "For a multi-dimensional array, accumulate is applied along only one\n", "axis (axis zero by default; see Examples below) so repeated use is\n", "necessary if one wants to accumulate over multiple axes.\n", "\n", "Parameters\n", "----------\n", "array : array_like\n", " The array to act on.\n", "axis : int, optional\n", " The axis along which to apply the accumulation; default is zero.\n", "dtype : data-type code, optional\n", " The data-type used to represent the intermediate results. Defaults\n", " to the data-type of the output array if such is provided, or the\n", " the data-type of the input array if no output array is provided.\n", "out : ndarray, optional\n", " A location into which the result is stored. If not provided a\n", " freshly-allocated array is returned.\n", "\n", "Returns\n", "-------\n", "r : ndarray\n", " The accumulated values. If `out` was supplied, `r` is a reference to\n", " `out`.\n", "\n", "Examples\n", "--------\n", "1-D array examples:\n", "\n", ">>> np.add.accumulate([2, 3, 5])\n", "array([ 2, 5, 10])\n", ">>> np.multiply.accumulate([2, 3, 5])\n", "array([ 2, 6, 30])\n", "\n", "2-D array examples:\n", "\n", ">>> I = np.eye(2)\n", ">>> I\n", "array([[ 1., 0.],\n", " [ 0., 1.]])\n", "\n", "Accumulate along axis 0 (rows), down columns:\n", "\n", ">>> np.add.accumulate(I, 0)\n", "array([[ 1., 0.],\n", " [ 1., 1.]])\n", ">>> np.add.accumulate(I) # no axis specified = axis zero\n", "array([[ 1., 0.],\n", " [ 1., 1.]])\n", "\n", "Accumulate along axis 1 (columns), through rows:\n", "\n", ">>> np.add.accumulate(I, 1)\n", "array([[ 1., 1.],\n", " [ 0., 1.]])\n" ] } ], "source": [ "np.info(np.add.accumulate)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`np.add.accumulate` 相当于 `cumsum` :" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": true }, "outputs": [], "source": [ "result_np = np.add.accumulate(y) * (x[1] - x[0]) - (x[1] - x[0]) / 2" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuczHX7x/HX5Uxko9JJ1F39UlTuTiLZKOfQ3VF0l3QT\nSSpRIcqhYpXSXblFUSJ0Ylch2ZwiwkaO1a7oQLLjnMP6/P74zmYas3bX2N2Znffz8ZjHfmfnMzOf\nx3fqcu31ub6fMeccIiISO4oU9ARERCR/KfCLiMQYBX4RkRijwC8iEmMU+EVEYowCv4hIjAkr8JtZ\nZTObbWbfmdlKM+uaxbhXzGy9maWYWc1w3lNERMJTLMznHwAecc4tN7OywDdmNtM5tzpzgJk1Bc5z\nzp1vZlcDrwO1wnxfERE5RmFl/M6535xzy/3Hu4DVwBlBw1oAY/xjFgFxZlYpnPcVEZFjd9xq/GZW\nFagJLAp66ExgY8D9TcBZx+t9RUQkd45L4PeXeSYDD/sz/yOGBN3XPhEiIgUk3Bo/ZlYc+AB41zn3\ncYghPwOVA+6f5f9d8OvoHwMRkVxyzgUn1tkKt6vHgFHAKufcsCyGTQH+7R9fC/A55zaHGuici8pb\n3759C3wOmn/Bz0Pzz6dbYiIuPd07Tk/Hde6MS0vzfp+WhqteHZeSgmvf3rulpuI6d2bb3JUsoBZj\nErbQ68F07uJdrq25ixOL96B48UPEldjFBeV/5Zp/7qXxWd9y57kLaX/7Dh6sNovHqn9Gzwd89Lxk\nGj0vmcbk1zZ775s5jwK6HatwM/46QFvgWzNb5v/dU8DZ/kA+wjk3zcyamtn3wG6gXZjvKSKxJikJ\n6tSBuDjvZ69e0KMHrFzp/WzeHMaNg8GDce+OY+OlzVny+hKWrCjB8nMWsKLSS/hGHuLCS5K54PVl\nnH/iFhoPqMvZo7rwUduTeOG57ZTs1d17rz59oP/LAcfvecdPXAj9P/COW18DDPTmMXCgN68oElbg\nd87NIwd/NTjnuoTzPiISg3IY7DOeH0JK7w+Zc+l/mXfDQOZdUwRXMZUrB63kypPTeODlS6nx3xZU\n+WoERU4qD93f9F7/vL0wux+zb72VksnTISHB+/2wYdkfz58PzZp5QT/zOJoU9J+JAX+yuGg1e/bs\ngp5CWDT/gqX5B0hMdC493TtOT3euc2fn0tK836elOVe9unMpKe6nu59yb/Te6G5loqtQ/oC7MO5X\n1/HWre6ds55wP87d5A5tS3eufXvvNmGC99zOnb3j9HTvlpjozX/q1L+Oo40/buY63poLo050PJmZ\ni5S5iEg+Cszsfb6/Z/bVq0Pz5rh3x7Gk/6d8FNeOxFG/8UvcRTQ+ZSkN7zmNBu+048zpow+PHzwY\nrrsOGjXyXj8zI/f5ojM7Pwozwx3D4q4Cv4gUrMxgn1kr37Dhr2C/8NkZTCh1Lx++t5cTqp7CzSd9\nSYtH/sFVL9xC0aQpfw/2mc8vhAE+Kwr8IhI9ArN88IJ19+5QuzarZ25ibJkHGD96D6XPOY3WJ3/O\nbU+dT7U+t0JiYswH+0DHGvi1O6eI5I+kJC9Aw+HF2g0bICmJHTvgjRV1qNX+IhrMepKMpSl8MnE/\nq8pcydOTalCt+Pde0B882HtulSqHF1bBC/4xFvTDoYxfRPJONvX7JQ168MYlr/NBUkkanLGa9o9X\n5MbXbqZY0ifK7HNApR4RiTwh6vf7m7ZiUttPGP5yBr8VPZOOv/Sl3S07Oe3NAV5QV7DPMQV+ESl4\nwbV78Mo5Dz7ItgGv8cYDy3n1xyZc9PuXPDTsPJqvSaDolf+EBQu8HvnAmr+CfbZU4xeRgnGU2j0+\nHz/1GUnXE9/mvJplWV/pWqbXH8znKafS8s2bKPrE43DffV7Q79Xr8OuoZp+nFPhFJDyZwd7n8wK2\n/6radRn/4L5rVlNz6jOUWr2MlTN/460f61HjhbawcePhxdrM5wUu1kqeUqlHRHLvKO2Y62Zt5Flf\nV6ZPO0iXThk8tG8oFYb2Uv0+D6jUIyL5JzDL9/tx64nc074odT7txYVb5/HDt3vou6E9FZ7ucrh0\nozbMiKCMX0RyJoss/7eLGzBgZCXG/3gVD7XbzSNftqL8pxO8IB/c1SPHlTJ+ETn+jrJwu2sX9P26\nGRc/2pASbh9rVmbQr/kSL+irdh/RlPGLSNZC9OFnNGvB6FuS6JtQlvqnr2Fghw1UWTdT7ZgFQH38\nInJ8HGXhdvaUHXRb3ZG4dV8ztMksrnjv0b9flauSTr5SqUdEjo8QC7dpf5TjlvbluW9xZ/qcM47k\nN3/gijN+OfwclXSiijJ+Eckyy//zyrok/O9EXvruRrq130X35OaUnvaBFm4jhDJ+ETl2IbL8Gann\nU+OB2nzze2W++foQfZp+4wV9LdxGPWX8IrHqKO2Zj7x6Lot+P5dX71tG010TtXAboZTxi0juBGX5\nhw7ByJXXUOPRG6hSegsr1xSn6bCG2kenEFLGLxJLssjy1/+jMR2GV2fPwRKMfDCFSzYmKcuPAsr4\nRSR7QVl+RgYkLGvANU/F0/KMJSxYXYFL+t6sLL+QU8YvUthlkeWvOacJ7V6+jFLFDjDq4RWcu366\nsvwoo4xfREILkeUPXVafa3vX4+6qc5m18jTO7XmbsvwYooxfJBb4e+5T73ySe27bA8Dbj6Qoy49y\nyvhF5LDAzdUAVz6O0ac+wVXXlaTlmUuY/d2pyvJjmAK/SGEUUN7ZuhVuabKHl5/fw+wnZ/JYzS8o\nWtQ/ThdhxSSVekQKixCLuDNbj6bdwg60LjaJAfPrU/ICbbVQmKjUIxLrArL8/fuhx9OlaDfnXsb4\nWjJk1uVe0Adl+aLALxLVAmv5/oD+fYfB1Pm/31nzyRqW3/wsDVJHwYgRf6v5q5Yf2xT4RaJZUKvm\nhA+KU3vKE/w77Vk+qTuUk1/tB1Wrehl+0CZsErtU4xeJNiFq+Xu7PUm3tIf5YtlJTLz+DWq2qAwL\nFqhVs5DTN3CJxIqgxdl1S3Zw240+LvbNY8Qdsyn3xhB9K1aM0OKuSKzIXJzt1YuJr27h2rrQ+eJk\nxr35J+XKuiPHaRFXgijjF4kGQeWd/fuh+71bSRq/nUlN3uKf73VXlh+DCizjN7PRZrbZzFZk8Xi8\nmW03s2X+W+9w31Mk5gQs4m7aBPHX/Ela4gq+eWEW/zzjt8PjlOVLDoSd8ZtZXWAXMNY5VyPE4/HA\no865Ftm8jjJ+kUAhFnGT24zkrnmdeaj0SHp+dTNFztEFWbGswDJ+59xcID2bYbmemEjMC8jynYMX\n/1uSO7/4D2N3tOTJGfW9oA/K8iXX8mNx1wG1zSzFzKaZ2UX58J4i0c8f0Hf3eIY2rXbx7ktbWNRi\nEDekvqkLsiQsxfLhPZYClZ1ze8ysCfAxcEGogf369fvrOD4+nvj4+HyYnkgECSrvpPniaDVnEJeu\nncj8O7+k9Osv/q2rR+Wd2JKcnExycnLYr3NcunrMrCowNVSNP8TYVOBy59y2oN+rxi8SUK//Ymkc\nd92RwZPFhtD1qbJYynJdkCV/U6AXcB0t8JtZJWCLc86Z2VXAROdc1RDjFPglNgVl+S7dx6tNpzHw\nu5a8V7Er9ZOfhipaxJUjFVjgN7PxQD3gZGAz0BcoDuCcG2FmDwKdgIPAHrwOn4UhXkeBX2JTQEDf\nXyaOzvfvY9HULUzxXcc5KZ/AJZf8fayyfPHTlg0i0cznY8sjz3HLymeouGU179R7k3LPPg5DhijD\nlyxpywaRaBL01Yjf/hTH1TP6U29JAh/WfYlyrwzUrpqSZxT4RQpCQI/+1KnQ4PoMBhbpw4DhcRQp\nVeLwOPXoSx5QqUckv4RYxB3aeAYvrW3KhxU7cPUXz2kRV3JFNX6RSBcQ0A+c4C3ifj1lM1N913J2\nSqIWcSXXVOMXiXT+sk36YwNocv1efp29hnk3vcDZqXN0Ja7kKwV+kbwUtIj747Y4as8eQPX5I/ik\n7lAt4kqBUOAXyUsBi7gLF0KdWgfp8udQhg0vRtHSWsSVgqEav0he8/mYfOdkOn19L2+f2JVmX/bU\nIq4cF1rcFYkUAd07zsHQoTBsyH6mbrmamiljtIgrx40Wd0Uihb+8c3Crjy5dYMyb+1kQ18wL+lrE\nlQigjF/keAjq0d/9s4/WtdPYW7oCk4vcQflPJ6i8I8edMn6RghSwiLtlC9RvWZaTiu0kae15lJ8w\nwgv6oEVciQgK/CLHgz+gr39wGLWvOkCjjE95O34MJVLXqbwjEUelHpFjFVTeWbQIWjU/wDNbH6RD\n270wfLj3mMo7kkdU6hHJbwHlncREaN40g5FlutFh+CVQsuThcSrvSIRRxi8SDp+PN29Oos/K2/n4\nxHu00ZrkK/Xxi+SHoB79/v3h7f/t47Ofa3BBymT16Eu+UqlHJD/4yzsZf/jo3Bk+nriPBWUbeUFf\ni7gSJZTxi+TSn7/5uKt2GjtO/Qcf7riREz99X+UdKRAq9YjklYDyjs8HLVvCGWV8jPnsVEqkLFF5\nRwqMSj0iecVf3vll9Xauuw4uO3cH4zZe5wV9lXckCinwi4QSuI9+XBzr2j1Hnav2c1fVBQxbXIci\nSVO9TF/76EsUUuAXCSWgR/+bbyD+prL0rv4JT0ytg703TlswSFRT4BcJxR/QZ90zliaNMnjt4tdo\nf/FCSE1VeUeingK/SKagr0mc/Hkcred2ZtIf9Wl1+iJISNDXJEqhoMAvkimgvDNiBHR98CAzKtxJ\nveG3aQsGKVTUzikSwKX7GNR4DqN+bsSMMjdz3szX1aMvEUt9/CK5FbS75qFD8FiXP/l8yh5m/Hwx\np6dMV4++RDT18YvkVkBp58ABaNdmP19/sIk51/fj9NSvtIgrhZYyfoltPh97e/bjjh8GcXDtD0yq\n/zonvDxI++hLVFCpRyQngso727dDiwa7OOObqYwZsY8St7f6e5BXeUcimEo9IjkR9N2419fZR/XU\nRMYtu5gSKYuPHK/yjhRCyvgl9vh8/NQ1gYZze3Pb7jE8+3VjrKo6dyT6qNQjkpWg8s7atdDw+v10\n+7Unj6S0U+eORC2VekSyElDeWboU4usepF+R/l7QV+eOxKCwAr+ZjTazzWa24ihjXjGz9WaWYmY1\nw3k/kRwL2l2TgQOZ02YEjevt4bXSj9Fu/v3aXVNiVrgZ/1tA46weNLOmwHnOufOBDsDrYb6fSM4E\nZPkASdOLceusTozfdRM3T22v3TUlpoUV+J1zc4H0owxpAYzxj10ExJlZpXDeUyRHMgN6r16Mf+V3\n2rfLYGrD4TRIHaXyjsS8vK7xnwlsDLi/CTgrj99TYlXQ7prExfH6SU/y+MP7+LzB81w99kHtrilC\n/izuBq84q3VH8kZAecc5GNQjnYTBjjl9v6B6pd8Pj1N5R2JcsTx+/Z+BygH3z/L/LqR+/fr9dRwf\nH098fHxezUsKI39Ad0/1oseefnw2wce8+aU5/cp/g6/F33v0Vd6RKJScnExycnLYrxN2H7+ZVQWm\nOudqhHisKdDFOdfUzGoBw5xztbJ4HfXxS+4F9ehnZEDH27ax4qP1fDqnLBXqXnx4rHr0pZApkD5+\nMxsPLAD+z8w2mtl9ZtbRzDoCOOemAT+a2ffACKBzOO8ncoSA8s6+fXBni92kzVjHrK9OoMKE17SI\nKxKCrtyV6OfzsbvHM/xr5bOUWfk14xefR6n/0xYMUvjpyl2JHUHdO+kujoaLB3DaVx8yKfkUL+iD\nFnFFsqDAL9EnoLyzeTPE197PlRs/5K1lNSk2Sj36ItlRqUeik89H2kNDuXFOb+7+8036LGqOdtiU\nWKPdOaVwC+reWbUKGtXfz+Obu9M15X7tsCkxSTV+KdwCyjuLF0P9ehkMKtbXC/ragkEkV5TxS/Tw\n+Zh979vcMacLb5Z7hBZzunubram8IzFKpR4pfILKOx9/DB3uO8jE9BuIT3lF5R2JeSr1SOETUN55\n+23o1OEgn1Zo4wV9lXdEjpkyfoksQVk+Ph8vNprOy+ubMP3ktlw4c7jKOyJ+yvilcAjaYfOpZ0vx\nvzV1mZtenQsnD9AXqIgcB8r4JfL4fGQ82ZvO259j6cytTGv4MqcM7AZDhijDFwmgxV2JXkHlnX37\n4O4WPrbOWMond4yn3BtDvMdU3hH5G5V6JHoFlHd27YLmN+wlY+Fipr24lnJlA5IBlXdEjgtl/BIZ\nfD62PjqIZov7ccnP03hjyZUUPVeLuCJHo1KPRJeg8s5PP0Gj+D9pmTqM55Y3xS5Vj75IdlTqkegS\nUN5ZtQquveYgHfa+wvMpTbH/qUdfJC8p45eC4/OxsP1IWn3ZjYTST9N23gPq0RfJBZV6JPIFlXem\nTYN72x7k7fQWNE15XlswiOSSSj0S+QLKO2PHwn33ZDDlpHu8oK8tGETyjTJ+yVtBWb5L95HQaCav\nft+Izyq2pdrn2oJB5Fgp45fIFJDlHzoEj/YqxZh11zA//SKqfaAtGEQKggK/5C1/QN/3RF/atNrF\nNx/+xNybBnNW6jyVd0QKiAK/HH9JSX8L6NstjqbLBrBv6gymX/88Jw1/FqpW9TJ8/18DIpJ/FPjl\n+Aso7/zyC1xXaz8Xfp/EpJd/pfQJAf/JqbwjUiC0uCt5w+djdadXaDLnCTpmvMYTC2/GqmoRV+R4\nUh+/FKyg7p25c+HWVgcYvO1+7kl5TD36InlAXT1SsALKO5MmwS03Z/BOuQe9oK9FXJGIooxfjl2I\nHv2XGk/nxXXNSKx4L5fNGqoefZE8pFKP5L+AgJ5RLo5unfYxe+LvJG2vQ5WUqSrviOQxlXok//m7\ncvb06MctTfawKvFH5t30AlVSv1R5RySCKfBL7gT16P/2ZxzxCwZRfuYkPo1/gbjh/dWjLxLhFPgl\ndwIWcb/7DmpdcZBmv7/F26/spESZYofHqUdfJGKpxi+55/Px+b/Hcte8zrx4Qh/toy9SQFTjl7wT\nVN7538Q42szvxKT0BrRNaq2N1kSijAK/ZM9f3sn4w0f37pDw/AHmVWxFvZThWsQViUIq9UhoQT36\nuzb5aFMnlR0lT+GDondQ4bP3VN4RKWAFVuoxs8ZmtsbM1ptZzxCPx5vZdjNb5r/1Dvc9JR8ELOL+\n9BPUaVyOU0v4mL7+XCq8/7rKOyJRLKzAb2ZFgVeBxsBFQGszqxZi6JfOuZr+24Bw3lPyiT+gL2w/\nklpXHuTespP533XjKJG6TuUdkSgXbsZ/FfC9cy7NOXcAmAC0DDEu13+KSAEIWsQd80kcLWY/wsgt\nLXnk/ERsaIJ69EUKgXAD/5nAxoD7m/y/C+SA2maWYmbTzOyiMN9T8krQIm7/vgdIrngLzYY3gZIl\nD49TeUckqhXLfshR5WQ1dilQ2Tm3x8yaAB8DF4T5vnK8BC7ixsXh6zGI1pf9xP4yO1lUujUVPxsX\nehFX5R2RqBVu4P8ZqBxwvzJe1v8X59zOgONPzew1M6vgnNsW/GL9+vX76zg+Pp74+PgwpyfZylzE\nHTiQVb/E0apFWZqc8AsJa/9J8ZRvQi/iKuCLFIjk5GSSk5PDfp2w2jnNrBiwFmgA/AJ8DbR2zq0O\nGFMJ2OKcc2Z2FTDROVc1xGupnTO/BLVq4vMx5bZ3uH/R/bxw2QTaXTAfeveGIUPUpikSwQqkndM5\ndxDoAkwHVgHvO+dWm1lHM+voH3YrsMLMlgPDgDvDeU85DgJaNQ8dgr7PleLBBXcxdWc87ap8AQla\nxBUpzHQBV6zy+fB1H0Db1GfZ8d0mJtZ/g9NuqA4LFniBP+CvAZV3RCKTvohFji6ovPPtt3DLTfto\n8tMIht61lOL/HeY9pitxRaKGNmmTowso77zzDjS4PoN+h/ryyvAiFC+t7ZRFYoky/sIsKMvft9nH\nI/HL+PyPmnxQ/j5qfP6S9tsRiWLK+OVIAVl+airUaVyWzbvLsvj3KtT4oJ/22xGJUQr8hU3gtgv+\ngD7ltne4+tK9tC05mck3jqB8aor22xGJYQr8hU1Alr9/PzzaqxQPfdWaT3bWp9v5SdpvR0QU+AuF\nEFl+asfnufbC3/khcTXLbu7PNaP+o/12RATQ4m7hELQ4+/7o3TzUOYNe+/rQtc027NXhatUUKYTU\nxx9rQmy7sPvhp+j6Yzfmfnsi468fyeUtztQFWSKFmAJ/rAnK3pfM3kmblruovfMzht8xn7JvJCjL\nFynk1M4Za/w1+ownezPo8XSaNnU8c/kU3hrlKFv2yHGq5YtIJmX80SK4tAP8uHwH99y+h2LrVzP2\n5o+pPPoZZfkiMUQZf2EU2K2T2aa5YQMuMYmRL+/hqlpGy5PmMGtkKpUr7D78PGX5InIUyvgjWXDm\nvmEDvzRqR4dTP+bXVemMrTeKi0c9qixfJEZpcbewCNGtQ/fuuGtq8857Rem+rA2d0gfSa0QVStze\n6u9BXh07IjFFpZ7CIuDK20ybtpXhpvtP5cWfb2fGjUN4JvUeSqQsPvK52nZBRHJAGX8kyCLLP1Sr\nNiPeKsHTS1vyULvdPPFlE0pM+1g7aooIoIw/uoXI8ldvqUj8f87jnXVX8+Vsx9PNvvGC/uDB3jgt\n4IrIMVLGX1CyyPL3XlGXgW9UZMTaevT910o6lRxN0ReH6MpbETmCMv5oEyLL/yz1Amp0qsO69FNI\nWXqILuOu8YJ+4DjV8UUkTMr481MWWX7aBQ155LXzWLHtLF65bzlNd03U/joiki1l/JEqi4uwSEpi\nzx54ZnFTruhZn8vjfmTluhI0HdbQC/rK8kUkjyjw57XAkk5cHPTogWvWnPe/q061CzJYtbMy37ww\ni95XfEapUv7naOFWRPKQSj15IYuSDrVrsyBxG91/7MSfKWt5uWESdd9/SFfeisgxUamnoB2lpAOw\n7veTuLX9idz5VVceqDiZJSOXU7fyhsPPV5YvIvlEGf/xEmJfHZo3Z9OLE3n2oS18tOGfPNphF90+\nv4nS0z7QRVgiEjZl/AUhxHfd0r07jB7Nb/3e4NHLvuDShqdSodgO1q46xJMNl3pBXxdhiUgBUuDP\nrWxKOpt9JenefhsXfdifjJSVrBw6g+drfUKFk5zXmVOlyt+DvTp2RCSfqdSTE4GLtZnlmR49YOVK\nqF4dmjfnp4SJDH7kV9774WruunkvTy6/gzOnj1ZJR0TyjEo9eSlESybNm0Plyqzo+S73nvslNRuf\nygm2h1XfHuTVuxd5QV8lHRGJQMr4s3KUlsxDS5Yys8ajvNR5Hd+e0oCHzp1Gx7t2UuHbZF1xKyL5\nRl/EcjwcraRTpw47Oj7O2IklebVqAqV2buXhro673m9JyWkfqaQjIvlOpZ5jldVi7fz5f5V0luy4\ngP9c/S1VPh7GlzcOYKR1ZNmSDNpd/q0X9FXSEZEoEpuBPwfBfkuZqgxrs5iaB77mtruKcU6Z31iV\ncpBJD8+n7uxnsSGDveeqS0dEokzhLvUElm4yjwGmT4c5c47ozNk1cjxT+ixm/P5bmDvnEC2aO+4t\nMpb4m8pRZOEC1e9FJKKoxp8pqzr9woUwc6Y3JiEBtm+H5s3ZMWI8Sf2+5iP7F9NnGHWuNVqX/IhW\n/S+nXIfWkJio+r2IRKQCq/GbWWMzW2Nm682sZxZjXvE/nmJmNY/5zQJLNJnHPh/063f4eNeukKUb\nqlX762V+SNnF8LaLaHziAs6qczbv/nkrjfYn8kPyJqb56nD3qHjKbdvgBX3V70WksHHOHfMNKAp8\nD1QFigPLgWpBY5oC0/zHVwMLs3gt59LTvVvfvqGPJ0xwrnNn59LSvOP27b1bWtrh4/R073716s6l\npDjXubPb+uVKN5l/uQfabHfnV93vKvGra9fyDzep0Ui3ff4Kb2xamnOJid7Pzp2913HO+5mY6ERE\nIo0XwnMfu4uF+e/GVcD3zrk0ADObALQEVgeMaQGM8f8js8jM4sysknNu8xGv1r2797NPn9DHCQlQ\nq5aXwY8bd7h0E1AiyvjDx5re77L49mksuPRT5p07lE0jHdfGv0OD9Yk8cNkqanz0L4rc3cbL6Feu\nPJzZZ5ZxMjP7Zs20WCsihU64gf9MYGPA/U14WX12Y84Cjgz8mQJr/YHHPh8MGQLjxuEuvZQtX29g\nXWpxVp7zLCvavsiKNSVYft4BTjv7ca5YsoLaPZvRafLd1FidQLGK5aH7DO911q5WsBeRmBVu4M/p\namzw4kPI5/UrW5ZDh2DfOdW47OWxVDuvDlvOuZfNQ8exeWtRNp0zlp/qD2LDDdv5vlwGxevu5vwT\nt3Dx7YOpsT6RW0/5mZpzm3BSpzv92fxy6JTgBfjrrvP+YgAvwAe2YSrYi0gUSE5OJjk5OezXCaur\nx8xqAf2cc439958EDjnnXggY8waQ7Jyb4L+/BqgXXOoxM1fUDuIwTiznKH8onfIl9nJKtZOp9FsK\nlUrt4KzG1Tl7YgJnD3+cf9QoQ8VBj3lP7tMH+vf3jm+80SsHBWbzar0UkUKoQNo5zawYsBZoAPwC\nfA20ds6tDhjTFOjinGvq/4dimHOuVojXcvvbdaCYZWBPBwTyrIL6dddBo0be74cNg27dvOPMAK9g\nLyKFXIH18ZtZE2AYXofPKOfcc2bWEcA5N8I/5lWgMbAbaOecWxridZxLT/fuBAZyBXURkZB0AZeI\nSIzRJm0iIpIjCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwi\nIjFGgV9EJMYo8IuIxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYo8IuI\nxBgFfhGRGKPALyISYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYo8IuIxBgFfhGRGKPALyIS\nYxT4RURijAK/iEiMUeAXEYkxCvwiIjFGgV9EJMYUO9YnmlkF4H2gCpAG3O6c84UYlwbsADKAA865\nq471PUVEJHzhZPxPADOdcxcAs/z3Q3FAvHOuZmEN+snJyQU9hbBo/gVL8y9Y0T7/YxFO4G8BjPEf\njwFaHWWshfE+ES/a/8PR/AuW5l+won3+xyKcwF/JObfZf7wZqJTFOAd8bmZLzOw/YbyfiIgcB0et\n8ZvZTOCyal8SAAAEA0lEQVS0EA/1CrzjnHNm5rJ4mTrOuV/N7BRgppmtcc7NPbbpiohIuMy5rOJ1\nNk80W4NXu//NzE4HZjvnLszmOX2BXc65oSEeO7aJiIjEMOdcrkvpx9zVA0wB7gFe8P/8OHiAmZUB\nijrndprZCUBD4JlQL3YskxcRkdwLJ+OvAEwEziagndPMzgBGOueamdm5wIf+pxQDxjnnngt/2iIi\ncqyOOfCLiEh0ytcrd82ssZmtMbP1ZtYzizGv+B9PMbOa+Tm/7GQ3fzOLN7PtZrbMf+tdEPMMxcxG\nm9lmM1txlDGRfO6POv9IPvcAZlbZzGab2XdmttLMumYxLiI/g5zMP1I/AzMrZWaLzGy5ma0ys5BV\nhwg+99nOP9fn3jmXLzegKPA9UBUoDiwHqgWNaQpM8x9fDSzMr/kdp/nHA1MKeq5ZzL8uUBNYkcXj\nEXvuczj/iD33/vmdBlzmPy4LrI2y//5zMv+I/QyAMv6fxYCFwLXRcu5zOP9cnfv8zPivAr53zqU5\n5w4AE4CWQWP+uijMObcIiDOzrK4PyG85mT9E6MVqzmuhTT/KkEg+9zmZP0TouQdwzv3mnFvuP94F\nrAbOCBoWsZ9BDucPEfoZOOf2+A9L4CVx24KGROy5hxzNH3Jx7vMz8J8JbAy4v8n/u+zGnJXH88qp\nnMzfAbX9fypOM7OL8m124Yvkc58TUXPuzawq3l8vi4IeiorP4Cjzj9jPwMyKmNlyvItNZzvnVgUN\niehzn4P55+rch9POmVs5XUUO/lcrUlafczKPpUBl59weM2uC1+J6Qd5O67iK1HOfE1Fx7s2sLDAZ\neNifOR8xJOh+RH0G2cw/Yj8D59wh4DIzKw9MN7N451xy0LCIPfc5mH+uzn1+Zvw/A5UD7lfG+1f1\naGPO8v8uEmQ7f+fczsw/yZxznwLF/W2v0SCSz322ouHcm1lx4APgXefcEde9EOGfQXbzj4bPwDm3\nHUgCrgh6KKLPfaas5p/bc5+fgX8JcL6ZVTWzEsAdeBeBBZoC/BvAzGoBPnd4P6CClu38zaySmZn/\n+Cq8dtlQtbhIFMnnPluRfu79cxsFrHLODctiWMR+BjmZf6R+BmZ2spnF+Y9LAzcCy4KGRfK5z3b+\nuT33+Vbqcc4dNLMuwHS8xYlRzrnVZtbR//gI59w0M2tqZt8Du4F2+TW/7ORk/sCtQCczOwjsAe4s\nsAkHMbPxQD3gZDPbCPTF606K+HMP2c+fCD73fnWAtsC3Zpb5P+1TeBdARsNnkO38idzP4HRgjJkV\nwUt233HOzYqW2EMO5k8uz70u4BIRiTH66kURkRijwC8iEmMU+EVEYowCv4hIjFHgFxGJMQr8IiIx\nRoFfRCTGKPCLiMSY/weoVZxsAST89wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = plt.plot(x, - np.cos(x) + np.cos(0), 'rx')\n", "p = plt.plot(x, result_np)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 速度比较" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "计算积分:$$\\int_0^x sin \\theta d\\theta$$" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import sympy\n", "from sympy.abc import x, theta\n", "sympy_x = x" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = np.linspace(0, 20 * np.pi, 1e+4)\n", "y = np.sin(x)\n", "sympy_y = vectorize(lambda x: sympy.integrate(sympy.sin(theta), (theta, 0, x)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`numpy` 方法:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slowest run took 4.32 times longer than the fastest. This could mean that an intermediate result is being cached \n", "10000 loops, best of 3: 56.2 µs per loop\n", "-2.34138044756e-17\n" ] } ], "source": [ "%timeit np.add.accumulate(y) * (x[1] - x[0])\n", "y0 = np.add.accumulate(y) * (x[1] - x[0])\n", "print y0[-1] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`quad` 方法:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10000 loops, best of 3: 40.5 µs per loop\n", "result = 3.43781337153e-15\n", "number of evaluations 21\n" ] } ], "source": [ "%timeit quad(np.sin, 0, 20 * np.pi)\n", "y2 = quad(np.sin, 0, 20 * np.pi, full_output=True)\n", "print \"result = \", y2[0]\n", "print \"number of evaluations\", y2[-1]['neval']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`trapz` 方法:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10000 loops, best of 3: 105 µs per loop\n", "-4.4408920985e-16\n" ] } ], "source": [ "%timeit trapz(y, x)\n", "y1 = trapz(y, x)\n", "print y1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`simps` 方法:" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1000 loops, best of 3: 801 µs per loop\n", "3.28428554968e-16\n" ] } ], "source": [ "%timeit simps(y, x)\n", "y3 = simps(y, x)\n", "print y3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`sympy` 积分方法:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100 loops, best of 3: 6.86 ms per loop\n", "0\n" ] } ], "source": [ "%timeit sympy_y(20 * np.pi)\n", "y4 = sympy_y(20 * np.pi)\n", "print y4" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }