{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import scipy as sp\n", "import scipy.sparse\n", "from scipy.sparse.linalg import spsolve\n", "from scipy.sparse import csc_matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Poisson Equation\n", "This notebook solves the Poisson equation\n", "$$ \\frac{\\partial^2 u }{\\partial x^2} = f(x) \\\\\n", "a \\le x \\le b \\\\\n", "f(a) = \\alpha \\; \\mathrm{and} \\; f(b) = \\beta\n", "$$\n", "by solving the matrix equation \n", "$$\n", "A \\vec{u} = \\vec{b}\n", "$$\n", "where\n", "$$\n", "A = \\frac{1}{h^2} \\left(\\begin{array}{c}-2 & 1 & ...& & \\\\ 1 & -2 & 1 & &\\\\ & 1 & & & \\\\ .& & & & \\\\ .& & & & 1 \\\\ & & & 1 & -2 \\end{array}\\right)\n", "$$\n", "is the finite difference approxiation of the second derivative operator, and\n", "$$\n", "\\vec{b} = \\left(\\begin{array}{c} f(a+h) - \\alpha/h^2 \\\\ f(a+2h) \\\\ f(a+3h) \\\\ . \\\\ . \\\\ . \\\\ f(a+Nh) - \\beta/h^2 \\end{array}\\right)\n", "$$\n", "In the following code,\n", "$$\n", "f(x) = e^{-100x^2}\n", "$$\n", "representing an almost point-mass centred at $x=0$. $u$ will then represent the gravitational potential, given than $\\alpha = \\beta$ to preserve the symmetry." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'u')" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAF3CAYAAABntJB9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl41OW9N/73JzOTTPaQFUgCIcugqKiIqMgmEKpd1NrTVrvZ1hZtiwvQc13tdc7V5/d4rv7O83uOgFTUSlu3elqrXaxt7YGwCYgiuIs0k4VAFsjKZJ3JbPfvj5nBFIFs8517lvfrunKZzPp2I3nnc9/3V5RSICIiIiIiSjRJugMQERERERHpwDJEREREREQJiWWIiIiIiIgSEssQERERERElJJYhIiIiIiJKSCxDRERERESUkFiGiIiIiIgoIbEMERERERFRQmIZIiIiIiKihMQyRERERERECcmsO8B45efnq7KyMt0xiIgS2ltvvdWllCrQnSMa8fsUEZF+Y/0+FXNlqKysDIcPH9Ydg4gooYnIcd0ZohW/TxER6TfW71NcJkdERERERAmJZYiIiIiIiBISyxARERERESUkliEiIiIiIkpILENERERERJSQWIaIiIiIiCghsQwREREREVFCYhkiIiIiIqKExDJEREREREQJybAyJCJPikiHiHx4nvtFRH4mIvUi8r6IzDMqCxERJa7JfD8SkTtFpC74cWfkUhMRUSQYORl6GsCNF7j/JgBVwY/VAB43MAsRESWupzGB70cikgvgfwG4BsACAP9LRKYYmpSIiCLKsDKklNoLoOcCD7kFwLMq4A0AOSIyzag8RNHC3t4Pv1/pjkGUMCbx/ehTAGqUUj1KqdMAanDhUjVpg8Ne/M+HJ418CyKiqFd7qh8ftvZG5L107hkqBtA84uuW4G2fICKrReSwiBzu7OyMSDgiI9R39ONTD+/FC4ebR38wEUXK+b4fRfz71NMHmnDPc29H7IcAIqJo9B9//QjffOpNDHt9hr9XTBygoJTaqpSar5SaX1BQoDsO0YTt/kcnlAJ2/qNDdxQiCqNwfZ/6+nUzkZ1qwaYaexjTERHFjjcau7G/vgv3LK1Aitlk+PvpLEOtAEpHfF0SvI0obu2tC/zG+PWGbnh8fs1piCjofN+PIv59Kstqweol5dj5jw68c+K0kW9FRBR1lFLYuN2OwswUfO3amRF5T51l6GUA3wie4nMtgF6lFBdKU9xyun04eKwHZXlpGBj24u3j/EGHKEqc7/vRNgCrRGRK8OCEVcHbDPXNhWXIS0/GRk6HiCjB7K/vwptNPVizvBJWi/FTIcDYo7V/C+B1ALNFpEVE7hKRe0TknuBDXgHQCKAewC8AfN+oLETR4OCxbri9fvzwU7NhSpIzUyIiMtZEvx8ppXoA/AeAQ8GPB4O3GSo9xYx7llZgX10XDjZ2G/12RERRQSmFh7bbUZyTii9fXTr6E8LEbNQLK6XuGOV+BeAHRr0/UbTZa+9CijkJKy8uwrwZOdhr78K/fkp3KqL4N5nvR0qpJwE8aUSuC/n6dTPxi32N2FBjx+9WXwsRiXQEIqKI2vWPDrzX7MD/ue2yiOwVComJAxSI4sHeuk4smJULq8WEJVUF+LCtF90Dw7pjEVEUslpMWLO8Em8e68Fr9ZwOEVF88/sVNmy3Y2ZeGr5wVUlE35tliCgCWh1O1HcMYKktcMrUElsBlAqsjSUiOpcvX12K6dlWPLS9FoHhFRFRfNp25BQ+OtmH+1dUwWKKbD1hGSKKgL32wP6gJcEydGlxNqakWfCqnfuGiOjcUswm3LeiCu82O7CLx/ETUZzy+RU21thRUZCOW64456XcDMUyRBQBe+2dmJplRVVhBgDAlCRYVFWAfXVd/I0vEZ3XF64qwYzcNGysscPv558VRBR//vJeG+o6BrC22gZTUuT3R7IMERnM6/Njf30Xltjy/2kT9JKqfHT2D+PoyX6N6YgomllMSXhgZRWOtPVh25FTuuMQEYWV1+fHwzvsuGhqJj596TQtGViGiAz2XosD/S7vmSVyIaGvecQ2EV3ILVcUo6IgHZt22OHjdIiI4sgf325FU/cQ1lXbkKRhKgSwDBEZLjT5mTdjyj/dXpRlRcmUVBw92acjFhHFCFOSYG21Dfb2AfzlvTbdcYiIwsLt9WPzzjrMLclG9ZwibTlYhogM1tHngghQmJnyifumZlnR0cfjtYnowj596TRcNDUTm3fWwevz645DRDRpvzvcjFaHE+uqbVqvpcYyRGSwjv5h5GekwHyOoyKLsqzo6HdpSEVEsSQpSbCu2oZjXYP44zutuuMQEU2Ky+PDll11mD9zypnLjujCMkRksPY+1zmnQgBQkJnCyRARjUn1nCLMLcnG5h11cHs5HSKi2PXfB0+gvW8Y61fN1joVAliGiAzX3jeMoizrOe8ryrKif9iLIbc3wqmIKNaIBKZDrQ4nXjjcrDsOEdGEDLm9eHxPPRZW5OG6ijzdcViGiIzW0T+MoqxzT4ZCEyNOh4hoLJbaCjB/5hRs2VUPl8enOw4R0bg9c+A4ugbcWL/KpjsKAJYhIkN5fH50Dw6jIPP8kyEgsJSOiGg0IoJ1q2w41efCfx88oTsOEdG49Ls8eGJvA5bNLsBVM3N1xwHAMkRkqK6BYSiF806GQre393MyRERjs7AiHwsr8vD4nnousSWimPKr/cfgGPJgXXV0TIUAliEiQ4WWvxWdZzJUGLy9g5MhIhqH9ats6Bpw49nXj+uOQkQ0Jo4hN3617xhWzSnC3JIc3XHOYBkiMlBo+VvheSZDWalmpJiT0MHJEBGNw1Uzc7FsdgF+/moD+l0e3XGIiEa1dW8jBtxerIuSvUIhLENEBgqVnPOdJiciKMqycs8QEY3bumobHEMePPVak+4oREQX1DUwjKcPNOEzl03DRVOzdMf5JyxDRAbq6HMhSYC89OTzPqaQ1xoiogmYW5KD6jlF+MW+RvQOcTpERNHr53sa4PL48MDK6JoKASxDRIZq7xtGXkYKzKbz/69WlGVFez8nQ0Q0fuuqbeh3ebF1X4PuKERE59Te58Kv3ziOW68sRmVhhu44n8AyRGSgjn7XeU+SCynM4mSIiCbm4mlZ+MzcaXjqtSZ0D/DPESKKPo/urofPr3D/iirdUc6JZYjIQO19w+c9SS6kMNOKgWEvBod5RC4Rjd/alTa4PD48sbdRdxQion/ScnoIv33zBL44vwQz89J1xzknliEiA3X0u857klxIaHLEE+WIaCIqCzNw65XFeOZAE4/pJ6KosmVXPQSCNcujcyoEsAwRGcbj86N70H3mWkLnE7qfJ8oR0UTdv6IKXr/CY3u4d4iIokNT1yBefKsFX7lmBopzUnXHOS+WISKDdA0MQ6nzH6sdwskQEU3WzLx0fGl+CX5z8ARaHU7dcYiI8LOddTAnCb6/rEJ3lAtiGSIySOhQhMLM0Q5QsAYfz8kQEU1caBnKll31mpMQUaKr7+jHS++24s6FZWd+zolWLENEBgktexttMpRlNSPFnMRlckQ0KcU5qbhjQSlePNyME91DuuMQUQLbtKMOqRYT7l5SrjvKqFiGiAzSHlz2NtrR2iKCoiwrl8kR0aT94IZKmJIEm3fW6Y5CRAnq6Mk+/O39k/jW9bOQl3Hhn4GiAcsQkUE6+1xIEozpD4KirBROhoho0gqzrPjGdTPxp3daUN8xoDsOESWgjTV2ZFrN+O7i6J8KASxDRIZp7xtGfkYKTEky6mMLM6288CoRhcU9SytgtZg4HSKiiHu/xYGaj9rx3cXlyE6z6I4zJixDRAZp73eNul8opDArhcvkiCgs8jJS8K3ry/CX99pw9GSf7jhElEA2bLdjSpoF37q+THeUMWMZIjJIR9/wqCfJhRRlWTEw7MXAsNfgVESUCFYvrkCm1YxNNXbdUYgoQRxu6sGr9k7cvbQCmdbYmAoBLENEhunod435OMlQaeLx2kQUDtlpFnx3cTm2f9SOD1p6dcchogSwYbsd+Rkp+MZ1M3VHGReWISIDeHx+dA24Rz1JLiS0nI5L5YgoXL51fRly0izYUFOrOwoRxbkD9V14vbEb319WgbRks+4448IyRGSAroHQBVfHNxniiXJEFC6ZVgvuWVqBPbWdeOt4j+44RBSnlFLYUGPH1CwrvnLNDN1xxo1liMgA7X1ju8ZQSGg5HU+UI6Jw+sZ1M5GfkYwN27l3iIiMscfeibeOn8aa5ZWwWky644wbyxCRAUJ7f8Z6mlyW1QyrJQkd/ZwMEVH4pCWb8f1llTjQ0I0DDV264xBRnFFKYeN2O0qmpOJL80t1x5kQliEiA7T3h5bJjW0yJCIozLSemSgREYXLV66ZgalZVmzcbodSSnccIooj2z9qxwetvbhvRRWSzbFZK2IzNVGU6+hzIUkC1/sYq6KsFO4ZIqKws1pMWLO8EoePn8ar9k7dcYgoTvj9Cptq7JiVn47brizWHWfCWIaIDNA1MIy8jBSYkmTMzynITEHnACdDRBR+X5pfipIpqdhYw+kQEYXH3z44iX+c6scDK6tgNsVupYjd5ERRzDHkwZS08V1wLCctGX1Oj0GJiCiRJZuTcN+KKrzf0ouaj9p1xyGiGOf1+bFphx22ogx8du503XEmhWWIyACOIQ9yUpPH9ZycVAscQx7+1paIDHHblcWYlZ+OjTV2+P38c4aIJu7P77ahsXMQa1faxrUKJhqxDBEZwOH0IHvckyELvH6FQbfPoFRElMjMpiTcv6IK/zjVj1c+PKk7DhHFKI/Pj80763DJ9Cx86pKpuuNMGssQkQF6h9zISR1nGQpOkhxDbiMiERHhc5dPR1VhBjbV2OHjdIiIJuD3b7XgRM8Q1q+yISnGp0IAyxCRIRxOD3LGORkKTZIcQ9w3RETGMCUJ1lbb0NA5iD+/26o7DhHFmGGvD4/srMMVpTm4YXah7jhhwTJEFGbDXh+G3D7kpI1/zxAA9PIQBSIy0I2XTMWcaVnYvLMOHp9fdxwiiiHPv9mMtl4XfrhqNkRifyoEsAwRhV1vcLKTPd5lcsHydJrL5IjIQElJgnXVNhzvHsIf3mrRHYeIYoTT7cOW3fVYMCsX11fm6Y4TNixDRGHmCE52xrtMLofL5IgoQlZcXIjLS3PwyK56DHt5aAsRje65N46js38Y66ttcTMVAliGiMIuVGbGe7R2NpfJEVGEiAh+uMqGVocTvzvUrDsOEUW5wWEvHn+1AYur8nFNefxMhQCWIaKwC50GN97JkNVigtWSxNPkiCgiFlXmY0FZLrbsqofLw+kQEZ3f0wea0DPoxrpqm+4oYccyRBRmoWVy490zBASmSVwmR0SRICJYv8qGjv5hPPfGcd1xiChK9To9eOLVBqy4qBBXzpiiO07YGVqGRORGEakVkXoR+dE57p8hIrtF5B0ReV9EPm1kHqJICB2gMN7JUOg5Di6TI6IIuaY8D4sq8/HYngYMDnt1xyGiKPSr/cfQ5/JibRxOhQADy5CImAA8CuAmAHMA3CEic8562L8DeEEpdSWA2wE8ZlQeokhxON0wJQkyUszjfm52quVMmSIiioR1q2zoGXTj6QNNuqMQUZQ5PejGk/uP4aZLp+LS4mzdcQxh5GRoAYB6pVSjUsoN4HkAt5z1GAUgK/h5NoA2A/MQRYRjyIOcVMuETloJTIa4Z4iIImfejClYflEhtu5tRJ+Lv4whoo89sbcRg+74nQoBxpahYgAjj6hpCd420v8D4Gsi0gLgFQD3GpiHKCIcTg+yJ7BEDuCeISLSY121Db1OD36175juKEQUJTr7h/HMgSbccvl02IoydccxjO4DFO4A8LRSqgTApwH8WkQ+kUlEVovIYRE53NnZGfGQROPRG5wMTURoz5BSKsypiIjO79LibNx4yVT8av8xnB7kdJqIgMf3NMDt8+P+lfE7FQKMLUOtAEpHfF0SvG2kuwC8AABKqdcBWAHkn/1CSqmtSqn5Sqn5BQUFBsUlCg+H042ctPFdYygkO80Ct9cPl8cf5lRERBe2ttqGQbcXT+xt1B2FiDQ72evEcweP4wvzijErP113HEMZWYYOAagSkVkikozAAQkvn/WYEwBWAICIXIxAGeLoh2KaY8gzoZPkAGBKsERx3xARRdrsqZm4+fLpeOZAEzr7h3XHISKNHt1dD6UU7l1epTuK4QwrQ0opL4A1ALYBOIrAqXFHRORBEbk5+LD1AL4rIu8B+C2AbyquD6IYF1gmN7HJUGh5HfcNEZEO96+ogtvnx+N7GnRHISJNmnuG8LtDzfjy1aUozU3THcdw4z/7dxyUUq8gcDDCyNt+MuLzjwBcb2QGokjy+PzoH/ZOeDIUOniBZYiIdCgvyMBtVxbjuYPHsXpJOaZmW3VHIqIIe2RXHUQEa26I/6kQoP8ABaK40uec+AVXAZyZKPVymRwRaXLfiir4/QpbdtfpjkJEEXasaxB/eLsVX7tmZsL8MoRliCiMHMEylD2J0+QAToaISJ/S3DR8+epS/O5QM1pOD+mOQ0QRtHmHHcmmJHxvWYXuKBHDMkQURqESM9HT5M6UISfLEFE4iciNIlIrIvUi8qNz3D9TRHaKyPsiskdESkbc939F5IiIHBWRn8lErqgcY9Ysr4SI4JGd9bqjEFGE2Nv78ef32nDnwjIUZKbojhMxLENEYRRa3jbR6wylWkxINiVxMkQURiJiAvAogJsAzAFwh4jMOethDwF4Vik1F8CDAP4z+NyFCOxtnQvgUgBXA1gaoejaTMtOxVevmYHfv92CY12DuuMQUQRsqrEjPdmMu5eU644SUSxDRGH08WRoYmVIRJCdZuGeIaLwWgCgXinVqJRyA3gewC1nPWYOgF3Bz3ePuF8hcNmHZAApACwA2g1PHAW+t6wCFpPgZzu5d4go3n3Y2ou/f3gK3140C1PSJ7a6JVaxDBGF0ZkyNMGjtQPPtXAyRBRexQCaR3zdErxtpPcA3Bb8/PMAMkUkL3hB8N0ATgY/timljhqcNyoUZlpx58IyvPRuK+ra+3XHISIDbaqxIzvVgrsWzdIdJeJYhojCyDHkhgiQaZ34qfU5aRacHuJkiCjCfghgqYi8g8AyuFYAPhGpBHAxgBIECtRyEVl89pNFZLWIHBaRw52d8XPt8LuXVCDNYsKmHXbdUYjIIO+cOI2d/+jA6iXlEz4AKpaxDBGFkcPpQXaqBUlJE99fnZ2azMkQUXi1Aigd8XVJ8LYzlFJtSqnblFJXAvi34G0OBKZEbyilBpRSAwD+DuC6s99AKbVVKTVfKTW/oKDAqL+PiMtNT8Zdi2bhlQ9O4Uhbr+44RGSAjTV25KYn45sLy3RH0YJliCiMHEOeCR+eEJKTZkEvT5MjCqdDAKpEZJaIJAO4HcDLIx8gIvkiEvqe+GMATwY/P4HAxMgsIhYEpkYJsUwu5K7F5ciymrGphtMhonhzsLEb++q68P1lFUhPmfiqlljGMkQURg6nB9kTPFY7hHuGiMJLKeUFsAbANgSKzAtKqSMi8qCI3Bx82DIAtSJiB1AE4KfB238PoAHABwjsK3pPKfWXSObXLTvVgtVLyrHjaAfeOXFadxwiChOlFDbU2FGYmYKvXTtTdxxtErMCEhmkd8g94WsMheSkWeD0+ODy+GC1mMKUjCixKaVeAfDKWbf9ZMTnv0eg+Jz9PB+Auw0PGOW+ef0sPPlaEzbW2PHru67RHYeIwuC1+m68eawHD95ySUL/vMHJEFEYOZyeCR+rHRKaLPVxqRwRRYmMFDPuWVqOfXVdePNYj+44RDRJSik8tL0W07Ot+PLVpaM/IY6xDBGFUVj2DAWf72AZIqIo8vVrA1elf2h7LZRSuuMQ0STsru3Au80O3LeiCinmxJ0KASxDRGHj8yv0ucKwZyg4WeK+ISKKJqnJJvxgWQXePNaD1+q7dcchoglSSmHDdjtm5KbhC1eV6I6jHcsQUZj0uzxQCmGYDAXKlIPXGiKiKHPHNTMwPduKDTWcDhHFqm1HTuFIWx8eWFkFi4lVgP8EiMIkNMmZ7J6hM5MhLpMjoiiTYjZhzfIqvHPCgd21HbrjENE4+fwKG2vsqChIxy1XFOuOExVYhojCJFReJn+AQuD5vVwmR0RR6IvzSzAjNw0btts5HSKKMX99vw329gGsrbbBNIkLxMcTliGiMAkta8tOndyeocwUM0xJAoeTy+SIKPpYTEm4b0UVjrT1YduRU7rjENEYeX1+PLyjDhdNzcSnL52mO07UYBkiCpPeME2GRATZvPAqEUWxW6+YjvKCdGysscPn53SIKBb86Z1WHOsaxLpqG5I4FTqDZYgoTM7sGZrkAQqh1+CeISKKVmZTEtautMHePoC/vt+mOw4RjcLt9WPzzjrMLclG9Zwi3XGiCssQUZiEylB2GMpQdpqFe4aIKKp95rJpuGhqJh7eUQevz687DhFdwItvNaPltBPrqm0Q4VRoJJYhojBxON3ItJphDsMxlVPSkrlniIiiWlKSYG21Dce6BvGnd1p1xyGi83B5fHhkZz3mz5yCpbYC3XGiDssQUZj0DnkmvV8oJId7hogoBqyaU4TLirOxeWcd3F5Oh4ii0W8OnsCpPhfWreJU6FxYhojCxOH0nLlg6mRxmRwRxQIRwbpVNrScduKFw8264xDRWYbcXjy2pwELK/KwsCJfd5yoxDJEFCaOIXcYJ0PJ6B/2wsN1+EQU5ZbZCnDVzCnYsqseLo9PdxwiGuHZ14+ja2AY61fZdEeJWixDRGHicHqQZQ1PGcpONQP4+LhuIqJoJSJYX23DqT4XfvvmCd1xiCio3+XBE682YNnsAlw1M1d3nKjFMkQUJn1OL7LCcJIcgDOv0+/yhuX1iIiMtLAyH9eV5+HR3Q1wujkdIooGT73WhNNDHqyvnq07SlRjGSIKA6UU+pyesByrDXx8PDcnQ0QUK9avsqFrYBjPvt6kOwpRwusd8uAX+xoDh5yUZOuOE9VYhojCYNjrh9vnR1ZwedtkhSZDfSxDRBQj5pflYqmtAD9/tQH9Lv7ZRaTTL/Y1YmDYi3XcKzQqliGiMAiVlnDtGQq9Th9/oCCiGLJ+lQ2nhzx4+rUm3VGIElb3wDCeeu1Y8MLIWbrjRD2WIaIwCC1n4zI5Ikpkc0tyUD2nCFv3NfLyAESaPLG3EU6PDw+s5FRoLFiGiMIgNMEJ3wEKgeV2fU4eoEBEsWVdtQ39Li9+ub9RdxSihNPR58IzB5pw65XFqCzM0B0nJrAMEYVBqLRkWcOzZyjVYoLFJFwmR0Qx5+JpWfjM3Gl4cv8x9Ay6dcchSiiP7WmAz69w/4oq3VFiBssQURiEe5mciCDLauEyOSKKSWtXVsHp8eGJVxt0RyFKGK0OJ35z8AS+OL8EM/PSdceJGSxDRGEQ7mVyodfiaXJEFIsqCzNx6xXFeOb1JnT0u3THIUoIW3bVAQDWLOdUaDxYhojCINynyQHBMsSLrhJRjLp/ZRU8PoXHdnM6RGS0492DePFwC+5YUIrinFTdcWIKyxBRGPQ6PUi1mJBsDt//UllWM5fJEVHMmpmXji9eVYLfHDyBNodTdxyiuLZ5Zx1MSYIf3FCpO0rMYRkiCoM+pzdsF1wNyUq1oJ9liIhi2L3BTdxbdtdrTkIUv+o7BvDSO634xnUzUZhl1R0n5rAMEYVBn8sT1iVyQOAwBp4mR0SxrDgnFbcvKMULh5pxontIdxyiuPTwDjusFhPuWVqhO0pMYhkiCoNepydsJ8mFhE6TU0qF9XWJiCLpBzdUwpQk2LyzTncUorhz9GQf/vr+SXzr+jLkZaTojhOTWIaIwqDP5QnrSXJA4MKrHp+Cy+MP6+sSEUVSUZYVX792Jv70TgsaOgd0xyGKK5tq7Mi0mrF6MadCE8UyRBQGfU5v2CdDodfjUjkiinX3LKuA1WLCphq77ihEceP9Fge2f9SO7ywqR3ZaeH8GSSQsQ0Rh0Ov0IMsa5gMUgnuQeKIcEcW6/IwUfHNhGf76/kn841Sf7jhEcWFjjR05aRZ8e1GZ7igxjWWIaJL8foV+A5bJnZkMsQwRURxYvaQcmSlmToeIwuCt4z3YU9uJu5dUIDPMBzglGpYhokkacHvhVwj/AQqpnAwRUfzISUvGXYtnYduRdnzQ0qs7DlFM27DdjvyMZNy5cKbuKDGPZYhokkKTm3AfrR1adsc9Q0QUL769aBZy0izYWFOrOwpRzDrQ0IUDDd343rJKpCWHd4l+ImIZIpqkPqcXAMJ+0dWPl8l5w/q6RES6ZFktWL2kHLtrO/HW8dO64xDFHKUUNm63Y2qWFV+9ZobuOHGBZYhokkLL2MK9ZyiTBygQURz65sIy5GckY8N2ToeIxutVeycOHz+NNcsrYbWYdMeJCyxDRJMUWsYW7mVyyeYkpFpMPECBiOJKWrIZ9yytwIGGbrze0K07DlHMUEphY40dJVNS8aX5pbrjxA2WIaJJCpWVcB+gEHpN7hkionjztWtnoigrBRtraqGU0h2HKCbUfNSO91t6cd/yKiSb+SN8uBj6T1JEbhSRWhGpF5EfnecxXxKRj0TkiIj8xsg8REYwaplc4DXNXCZHRHHHajFhzQ2VONR0GnvrunTHIYp6fn9gKlSWl4bb5hXrjhNXDCtDImIC8CiAmwDMAXCHiMw56zFVAH4M4Hql1CUAHjAqD5FR+lxeiACZKeE/0SU71cIDFIgoLn3p6lIU56Ri43ZOh4hG88qHJ/GPU/14YKUNZhOnQuFk5D/NBQDqlVKNSik3gOcB3HLWY74L4FGl1GkAUEp1GJiHyBB9Tg8yUsxISpKwv3aWlcvkiCg+pZhNuH9FFd5r6cWOo/z2T3Q+Pr/Cpho7qgoz8LnLp+uOE3eMLEPFAJpHfN0SvG0kGwCbiLwmIm+IyI0G5iEyRJ/TY8h+ISCw9I7L5IgoXt02rxhleWnYWGOH38/pENG5/PndVjR0DmJttQ0mA37xmuh0z9nMAKoALANwB4BfiEjO2Q8SkdUiclhEDnd2dkY4ItGF9bk8YT9JLiSwTI5liIjik9mUhAdW2nD0ZB/+/uEp3XGIoo7H58fmnXWYMy0LN14yVXecuGRkGWoFMPLcv5KA6XWeAAAgAElEQVTgbSO1AHhZKeVRSh0DYEegHP0TpdRWpdR8pdT8goICwwITTUSf0xv2C66GZFnN6B/28jemRBS3Pnf5dFQVZmDTDjt8/LOO6J/84a0WHO8ewrpqmyHL8cnYMnQIQJWIzBKRZAC3A3j5rMe8hMBUCCKSj8CyuUYDMxGFXa/By+SUAvqHeYgCEcUnU5JgbbUN9R0DePm9s39nSpS4hr0+/GxnHS4vzcGKiwt1x4lbhpUhpZQXwBoA2wAcBfCCUuqIiDwoIjcHH7YNQLeIfARgN4B/VUrxCmwUU4xcJhc6rptL5Ygont14yVRcPC0Lm3fUwePz645DFBV+d6gZbb0u/HCVDSKcChnF0D1DSqlXlFI2pVSFUuqnwdt+opR6Ofi5UkqtU0rNUUpdppR63sg8REboc3oMucYQgDMliyfKEVE8S0oSrK+2oal7CH98u0V3HCLtXB4ftuyqx4KyXCyqzNcdJ67pPkCBKKZ5fH4Mun2GLZMLvS5PlCOieLfi4kJcXpqDn+2sx7DXpzsOkVbPvXEcHf3DWM+pkOFYhogmod8V2MuTZTXoAIXgwQy88CoRxTuRwHSo1eHE7w41j/4Eojg1OOzFY3sasKgyH9eU5+mOE/dYhogmIbSXJzuNy+SIiCZrcVU+FpTlYsuuerg8nA5RYnr6QBN6Bt1Yt8qmO0pCYBkimoTQ8jXDrjOUxgMUiChxiAjWrbKho38Yz71xXHccoojrdXrwxKsNWH5RIebNmKI7TkJgGSKahNDExqgDFDKSzRBhGSKixHFteR4WVebj8T0NGORlBSjB/Gr/MfS5vFhXzalQpLAMEU1CaDJk1AEKSUmCzBQzD1AgooSybpUN3YNuPPN6k+4oRBFzetCNJ/cfw02XTsWlxdm64yQMliGiSQgdbGDUMjkgsFSuz8XfjhJR4pg3YwqWX1SIJ15t5J5JShhP7G3EoNuLtZwKRRTLENEkfLxMzpjT5IBA0eIyOSJKNOuqbeh1evDk/mO6oxAZrrN/GM8caMLNl0+HrShTd5yEwjJENAm9Tg8sJkGqxWTYe2RZLVwmR0QJ59LibNx4yVT8ct8xOIbcuuMQGerxPQ0Y9vpw/4oq3VESDssQ0ST0OT3IsloMvSBadqqFy0SIKCGtrbZh0O3F1r2NuqMQGeZUrwvPHTyOL8wrQXlBhu44CYdliGgS+lxew06SC8lKNfOiq0SUkGZPzcTn5k7HU681oWtgWHccIkNs2V0Hv1/hPk6FtGAZIpqEXqfH8DKUncplckSUuB5YWYVhrw+P72nQHYUo7Jp7hvC7Q8348tWlKM1N0x0nIbEMEU1CYJmccYcnAIE9Q06PD26v39D3ISKKRuUFGbhtXgmee+M4TvW6dMchCqtHdtVBRLBmeaXuKAmLZYhoEvpcxk+GQq/fz31DRJSg7l9RBZ9f4dHd9bqjEIXNsa5B/OHtVnz1mhmYlp2qO07CYhkimoQ+p8ewC66GhF6fS+WIJk5EbhSRWhGpF5EfneP+mSKyU0TeF5E9IlIy4r4ZIrJdRI6KyEciUhbJ7ASU5qbhS1eX4vlDJ9DcM6Q7DlFYbN5hh8Uk+N6yCt1REhrLENEEKaXQ5/QaesFV4ONrGPHCq0QTIyImAI8CuAnAHAB3iMicsx72EIBnlVJzATwI4D9H3PcsgP9SSl0MYAGADuNT09nuXV4JgWDLLk6HKPbVtffjz++14c7rylCYadUdJ6GxDBFNkMvjh9vnN/SCqwDOlC1OhogmbAGAeqVUo1LKDeB5ALec9Zg5AHYFP98duj9YmsxKqRoAUEoNKKU4mtBgWnYqvnLNDPz+7RY0dQ3qjkM0KZt22JFmMeHupZwK6cYyRDRBoXKSk5ps6PvkpLEMEU1SMYDmEV+3BG8b6T0AtwU//zyATBHJA2AD4BCRP4rIOyLyX8FJE2nw/RsqYDEJNu+s0x2FaMKOtPXilQ9O4a5Fs5CbbuzPEDS6MZUhESkUkc+LyA9E5NsiskBEWKQooYXKidF7hrK4Z4goEn4IYKmIvANgKYBWAD4AZgCLg/dfDaAcwDfPfrKIrBaRwyJyuLOzM2KhE01hphV3XleGl95tRV17v+44RBOyqcaOLKsZdy0u1x2FMEoZEpEbRGQbgL8hsNZ6GgJLCf4dwAci8r9FJMv4mETRxzHkBvDx5MYoZw5QCL4fEY1bK4DSEV+XBG87QynVppS6TSl1JYB/C97mQGCK9G5wiZ0XwEsA5p39BkqprUqp+Uqp+QUFBUb9fRCAu5dWIM1iwsM7OB2i2PNuswM7jnZg9ZJyw3+ZSmMz2maHTwP4rlLqxNl3iIgZwGcBVAP4gwHZiKJapCZDKWYTUi0mToaIJu4QgCoRmYVACbodwFdGPkBE8gH0KKX8AH4M4MkRz80RkQKlVCeA5QAORyw5fUJuejK+vWgWHtlVj++39eKS6dm6IxGN2YbttchNT8Y3r5+lOwoFXXAypJT613MVoeB9XqXUS0opFiFKSI4IlSEgMH1yDLEMEU1EcKKzBsA2AEcBvKCUOiIiD4rIzcGHLQNQKyJ2AEUAfhp8rg+BJXI7ReQDAALgFxH+W6CzfGdxObKsZmyq4XSIYsebx3qwr64L9ywtR0aKsYcv0diN6d+EiPwawBqlVG/w6zIAv1JKrTAuGlF06wuVIYOXyQGBwsXJENHEKaVeAfDKWbf9ZMTnvwfw+/M8twbAXEMD0rhkp1rw3cXl2FBjx7vNDlxRmqM7EtEFKaWwYXstCjJT8PVry3THoRHGWkv3AzgoIusQOIHnXwGsNywVUQxwDHlgShJkRuC3O1mpljOTKKJEJiI/OdftSqkHI52F9PrWoll48rVj2Fhjx7PfXqA7DtEFHWjoxsFjPfjfN1+C1GQeSBlNxnQinFLqCQDfAfBnBC5Gt0Qp9RcjgxFFu16nB1lWM0TE8PfKSbWcmUQRJbjBER8+BA73KdMZiPTISDHjnqUV2GvvxKGmHt1xiM5LKYWHttdierYVty8oHf0JFFFjPVr76whsJv0GgKcBvCIilxuYiyjqOZwe5KRF5voA2ancM0QEAEqpDSM+forAXh+eT5ugvnFdGfIzUvDQtloopXTHITqn3bUdeOeEA/euqEKKmVOhaDPWawV9AcAipdRvlVI/BnAPAqWIKGH1Oj1nrgFktJw07hkiOo80BI7KpgSUmmzCmhsqcPBYDw40dOuOQ/QJSilsrLFjRm4a/uUq/lEVjca6TO5WpVTHiK/fBHCNYamIYkDvkBs5ESpD2akWOD0+DHt9EXk/omglIh+IyPvBjyMAagE8rDsX6XPHNTMwLduKh7ZzOkTRZ9uRU/iwtQ/3raiCxTTWGQRF0mgXXf13Eck9131KKbeILBeRzxoTjSi69To9EbtgWnZwOR6nQ0T4LIDPBT9WAZiulNqiNxLplGI24d7lVXjnhAN7ajt1xyE6w+cPTIUqCtLx+SuLdceh8xjtGKwPAPxFRFwA3gbQCcAKoArAFQB2APh/DU1IFKUCe4YiNxkCgN4hDwozrRF5T6JopJQ6rjsDRZ8vzi/B46/WY0NNLZbNLojIwTZEo/nr+22wtw/gkTuuhCmJ/01Gq9Hmdf+ilLoegQvVHQFgAtAH4DkAC5RSa4NX5CZKKH6/Ql8EJ0Oh5XicDBERfZLFlIT7V9jwYWsfth1p1x2HCF6fH5t31OGiqZn4zGXTdMehCxitDF0lItMBfBXAywCeAPAsgEMAUg3ORhS1+oe98CtEbplc8H14ohwR0bndesV0lOenY1ONHX4/9w6RXn96pxWNXYNYW21DEqdCUW20MvRzADsBXATg8IiPt4J/JUpIoWv+RGwylMbJEBHRhZhNSXig2oba9n789YOTuuNQAnN7/fjZrjpcVpyNVXOKdMehUVywDCmlfqaUuhjAk0qp8hEfs5RSvK4DJazQhCaS1xkCAvuUiIjo3D572TTMLsrEwzV2eH1+3XEoQb34VjOae5xYt8rG/WsxYKxHa3/P6CBEsaQ3wpOhTKsFIpwMERFdSFKSYG21DY1dg3jp3TbdcSgBuTw+bNlVj6tmTsEyW4HuODQGPPCcaAIcTjcAROw0OVOSIDPFfGZ5HhERndunLinCpcVZ2LzTDreX0yGKrN++eQIne11YX82pUKxgGSKagEhPhoDAkjzHkDti70dEFItEBOurZ6O5x4kX32rWHYcSiNPtw6O7G3BdeR4WVubrjkNjxDJENAGhPUORLEPZqRYukyMiGoNlswswb0YOtuyqh8vj0x2HEsSzrzeha2AY61fZdEehcWAZIpqAPqcHKeYkWC2miL1nTpqFBygQEY2BiGD9qtk42evCb988oTsOJYCBYS9+/moDltoKML8sV3ccGgeWIaIJcAxF7oKrIVmcDBERjdnCijxcW56LR3c3wOnmdIiM9dT+Yzg95MG6ak6FYg3LENEE9Do9ETs8ISQ71YJeXnSViGhMQtOhroFhPPt6k+44FMd6hzzYuq8RKy8uwuWlObrj0DixDBFNgMPpjvhkKCc4GVKKV1YnIhqLq8tyscRWgJ+/2oCBYa/uOBSnfrGvEf0uL/cKxSiWIaIJ6HV6kZ0amQuuhmSnWuD1KwxyuQcR0Zitq7bh9JAHT+0/pjsKxaHugWE89doxfGbuNFw8LUt3HJoAliGiCegd0jAZCi7L474hIqKxu6I0BysvLsLWfY1cakxh98TeRjg9PqxdWaU7Ck0QyxDRBOjaMwSA1xoiIhqnddU29Lu8+OX+Rt1RKI509Lnw7OtNuPWKYlQWZuqOQxPEMkQ0Th6fH4NuX8QnQ6FleZwMERGNz5zpWfjMZdPw5P5j6B4Y1h2H4sRjexrg8Snct4JToVjGMkQ0TqEyomsyxGUeRETjt7a6CkMeH7bu5XSIJq/N4cRvDp7Av8wrQVl+uu44NAksQ0Tj5AiWEe4ZIiKKHZWFmbj1imI883oTOvpduuNQjHtkVz0UFO5dUak7Ck0SyxDROIXKSOSXyQX3DLEMERFNyP0rquDxKTy2u0F3FIphJ7qH8OLhZtyxYAZKpqTpjkOTxDJENE69zsABBpEuQ2nJJlhMwskQEdEEleWn41/mleA3B0+gzeHUHYdi1OaddTAlCX5wA6dC8cDQMiQiN4pIrYjUi8iPLvC4L4iIEpH5RuYhCoeP9wxF9jpDIoLsVMuZZXpERDR+966ohILClt31uqNQDGroHMCf3mnB16+diaIsq+44FAaGlSERMQF4FMBNAOYAuENE5pzjcZkA7gdw0KgsROGka89Q6D37OBkiIpqwkilpuP3qGXjhUDOae4Z0x6EYs3lHHawWE+5ZVqE7CoWJkZOhBQDqlVKNSik3gOcB3HKOx/0HgP8PAHczUkwITYayrOaIv3d2qgUOJ68zREQ0GWuWV8KUJNi8s053FIoh/zjVh7+834ZvLixDfkaK7jgUJkaWoWIAzSO+bgnedoaIzANQqpT6m4E5iMLKMeRBZooZZlPkt9zlpCVzzxAR0SQVZVnxtWtn4o9vt6Cxc0B3HIoRm2rsyEg2Y/WSct1RKIy0HaAgIkkANgJYP4bHrhaRwyJyuLOz0/hwRBfQ5/QgO8LXGArhniEiovD43rIKpJhNeHgHp0M0ug9aerHtSDu+s7g84nuGyVhGlqFWAKUjvi4J3haSCeBSAHtEpAnAtQBePtchCkqprUqp+Uqp+QUFBQZGJhqdw+nRsl8ICJQhToaIiCYvPyMF37q+DH95vw21p/p1x6Eot7GmFjlpFnx7UZnuKBRmRpahQwCqRGSWiCQDuB3Ay6E7lVK9Sql8pVSZUqoMwBsAblZKHTYwE9Gk9To9Zy6AGmnZqRb0u7zw+ZWW9yciiierl5QjI9mMTTV23VEoir11/DR213bi7iUVyLTq+f5PxjGsDCmlvADWANgG4CiAF5RSR0TkQRG52aj3JTKaY8itdTIEgCfKERGFQU5aMu5aPAv/c+QUPmzt1R2HotTGmlrkZyTjzoUzdUchAxi6Z0gp9YpSyqaUqlBK/TR420+UUi+f47HLOBWiWNDr9CI7Vc964dBEikvliIjC49uLZiEnzYIN22t1R6Eo9HpDN16r78Y9SyuQlhz5U2TJeNoOUCCKRUop9Dr1T4YcLENERGGRZbVg9ZJy7K7txFvHT+uOQ1FEKYWNNbUoykrB167lVChesQwRjYPT44PHp7TtGeJkiIgo/O68rgx56cncO0T/ZG9dFw41ncaa5VWwWky645BBWIaIxiF0rLX2ydAQL7xKRBQu6SlmfG9ZBfbXd+H1hm7dcSgKKKWwcXstinNS8eX5paM/gWIWyxDROIQmMjnaylBgrxIPUCAiCq+vXTsTRVkp2FhTC6V4Ymei23G0A++19OL+FVVINvPH5XjGf7tE4xA9kyGWISKicLJaTFhzQyUONZ3Gvrou3XFII79fYWONHWV5abhtXrHuOGQwliGicTgdXJ42JV3PaXLJ5iRkpJjRw2VyRERh96WrS1Gck4oN2zkdSmR///AUjp7swwMrbTCb+KNyvOO/YaJx6B4MlJA8TWUIAHLTk9EzyDJERBRuKWYT7ltRifdaerHzaIfuOKSBz6+waYcdVYUZ+Nzl03XHoQhgGSIah54BvZMhgGWIiMhIt80rwcy8NGyoscPv53Qo0bz8XivqOwawttoGU5LojkMRwDJENA7dg8PITrXAonFsnp+RjK4BliEiIiNYTEl4YGUVjp7sw/8cOaU7DkWQx+fHwzvqMGdaFm68ZKruOBQhLENE49A96Na6RA4ITYaGtWYgIopnN19ejMrCDGysscPH6VDC+OPbLTjePYR11TYkcSqUMFiGiMahZ8CNXO1lKAU9g25u7iUiMogpSbB2pQ31HQP4y3ttuuNQBAx7ffjZznpcXpqDFRcX6o5DEcQyRDQO3YPD2stQfkYyPD6FPpdXaw4ionh206VTcfG0LGzaYYfH59cdhwz2u0PNaHU4sb7aBhFOhRIJyxDROPQMupGXkaI1Q6iM8RAFIiLjJCUJ1lfbcLx7CH96u1V3HDKQy+PDll31uLpsChZX5euOQxHGMkQ0Rn6/wukhT1TsGQLAfUNERAZbcXEhLi/NweaddXB7OR2KV8+9cRwd/cNYv2o2p0IJiGWIaIx6nR74/Er7Mrm89MBkqpsnyhERGUokMB1qdTjxu8PNuuOQAQaHvXh8TwMWVebj2vI83XFIA5YhojHqDk5i8jI0l6Hg+3dzmRwRkeEWV+Xj6rIp2LKrDi6PT3ccCrOnDzShe9CNdatsuqOQJixDRGMUmsSEJjO6cM8QEVHkiAjWr5qN9r5h/PfBE7rjUBj1uTzYurcRyy8qxLwZU3THIU1YhojGKFQ+dC+Ts1pMSE82cZkcEVGEXFueh+sr8/D4nnoMuXmSZ7z41b5j6HV6sK6aU6FExjJENEZdwTKke5lcIEPKmWV7RERkvHXVs9E14MbTB5p0R6EwcAy58eT+Y/jUJUW4tDhbdxzSiGWIaIx6gpOYKWn6y1BuejKXyRERRdBVM6fghtkF2Lq3Ef0uj+44NElb9zZiwO3FuurZuqOQZixDRGPUMziMLKsZyWb9/9vkpSdzmRwRUYStq54Nx5AHT+5v0h2FJqFrYBhPvdaEz82djtlTM3XHIc30/1RHFCO6o+CCqyGcDBERRd5lJdn41CVF+OW+RjiG+GdwrPr5ngYMe314YGWV7igUBViGiMaoe8Ct/fCEkNCeIaWU7ihERAllbbUNA24vtu5t1B2FJqC9z4Vfv3Ect80rQXlBhu44FAVYhojGqGfQjbxoKUPpyfD4FPqHeaoREVEkXTQ1C5+dOx1PH2hC1wAPsok1j+6uh8+vcP8KToUogGWIaIwCy+SiowydudYQ9w0REUXcAyur4PL48PM9Dbqj0Di0nB7Cb988gS9fXYrS3DTdcShKsAwRjYHfr3B6KJqWyQVy8HhtIqLIqyjIwOevLMGv3ziO9j6X7jg0Ro/srIeIYM3ySt1RKIqwDBGNQa/TA59fITc9Og5QyAvm4IlyRGMjIjeKSK2I1IvIj85x/0wR2Ski74vIHhEpOev+LBFpEZEtkUtN0ez+FVXw+RUe3V2vOwqNQVPXIH7/dgu+smAGpmWn6o5DUYRliGgMuoMnt+VHyzK5YA6eKEc0OhExAXgUwE0A5gC4Q0TmnPWwhwA8q5SaC+BBAP951v3/AWCv0VkpdszIS8MX55fit2+eQMvpId1xaBSbd9bBYhJ8/4YK3VEoyrAMEY1BqHREzTK59NAyOZYhojFYAKBeKdWolHIDeB7ALWc9Zg6AXcHPd4+8X0SuAlAEYHsEslIMuXd5JQSCLbs4HYpmde39eOndVtx5XRkKM62641CUYRkiGoPu4IlB0VKGrBYT0pNNXCZHNDbFAJpHfN0SvG2k9wDcFvz88wAyRSRPRJIAbADwQ8NTUsyZnpOKr1wzAy++1YKmrkHdceg8Ht5RhzSLCXcv5VSIPolliGgMQhOYvCjZMwQElsr18AAFonD5IYClIvIOgKUAWgH4AHwfwCtKqZYLPVlEVovIYRE53NnZaXxaihrfX1YBi0mweWed7ih0Dh+19eFvH5zEtxfNippfaFJ0YRkiGoNoWyYHALnpKVwmRzQ2rQBKR3xdErztDKVUm1LqNqXUlQD+LXibA8B1ANaISBMC+4q+ISL/5+w3UEptVUrNV0rNLygoMOhvg6JRYZYVd15XhpfebUV9R7/uOHSWjTV2ZFrN+M6ict1RKEqxDBGNQffAMDKtZiSbo+d/mfz0ZC6TIxqbQwCqRGSWiCQDuB3AyyMfICL5wSVxAPBjAE8CgFLqq0qpGUqpMgSmR88qpT5xGh0ltruXViDNYsKmHZwORZN3mx3YcbQdqxeXIzvNojsORano+cmOKIp1D7rPHFoQLXLTk3maHNEYKKW8ANYA2AbgKIAXlFJHRORBEbk5+LBlAGpFxI7AYQk/1RKWYlJuejK+vWgW/vb+SXzU1qc7DgVtrLFjSpoF31o0S3cUimIsQ0Rj0DMYPRdcDQnsGXJDKaU7ClHUU0q9opSyKaUqlFI/Dd72E6XUy8HPf6+Uqgo+5jtKqU9syFNKPa2UWhPp7BQbvrOoHJlWMzbW2HVHIQCHmnqw196J7y2rQEaKWXccimIsQ0Rj0DPoRl5G9ByeAAD56Slw+/wYGPbqjkJElPCy0yxYvbgcO462471mh+44CU0phYe21aIgMwVfv7ZMdxyKcixDRGPQNRCdy+QAcN8QEVGU+NaiWZiSZsEGToe0OtDQjYPHevCDZRVITTbpjkNRjmWIaBR+v8LpoehcJgfwwqtERNEiI8WMe5ZWYK+9E4eaenTHSUhKKWzYXotp2VbcvmCG7jgUA1iGiEbR5/LA51dRV4ZCkyoeokBEFD2+cV0Z8jNSsGF7re4oCWlPbSfePuHAvcurYLVwKkSjYxkiGkVXcBlafpTtGQrtYeoe4IVXiYiiRWqyCT+4oQJvNPbgQH2X7jgJRSmFDTW1mJGbhi/OL9Edh2IEyxDRKKLxgqvAx5MhLpMjIooudyyYgWnZVjy0vZYnfkbQtiPt+LC1D/etqILFxB9xaWz4XwrRKHoGA5OXaCtDVosJackmLpMjIooyVosJa5ZX4u0TDuyxd+qOkxD8foVNNXaUF6Tj1ium645DMYRliGgUoclLtC2TA4C8DF54lYgoGn3xqlKUTEnFxu12Toci4K8fnERtez8eWGmDmVMhGgf+10I0itDR1VPSLZqTfFJuegq6uGeIiCjqJJuTcP+KKnzQ2ottR9p1x4lrXp8fD++wY3ZRJj572TTdcSjGsAwRjeJkrwtT0ixIMUffqTRFmSk41evSHYOIiM7h81cWozw/HZtq7PD7OR0yykvvtqGxcxBrq21IShLdcSjGsAwRjeJkrxPTslN1xzin6TmpOMkyREQUlcymJNy/sgq17f342wcndceJSx6fH5t32nFpcRY+dUmR7jgUg1iGiEZx0uHC9JzoLEPTsq0YGPaiz+XRHYWIiM7hc3OnY3ZRJjbtsMPr8+uOE3dePNyC5h4n1lfPhginQjR+LENEo2jrdWJ6jlV3jHOaFixpJx2cDhERRaOkJMHa6io0dg7iz++26Y4TV1weHx7ZVYd5M3KwbHaB7jgUo1iGiC5gYNiLfpc3epfJZQdKWluvU3MSIiI6n09dMhWXTM/C5p118HA6FDbPv3kCJ3tdWL+KUyGaOJYhogs46QiUDE6GiIhookQEP1w1Gyd6hvDi4RbdceKC0+3Do3sacG15LhZW5OmOQzGMZYjoAtqChxNE62SoKDMFSRI45IGIiKLXstkFuHJGDrbsqsOw16c7Tsx79vUmdPYPcypEk2ZoGRKRG0WkVkTqReRH57h/nYh8JCLvi8hOEZlpZB6i8QpNhqZlR+dkyGxKQmGmFW2cDBERRbXQdKit14Xn32zWHSemDQx78fNXG7DEVoCry3J1x6EYZ1gZEhETgEcB3ARgDoA7RGTOWQ97B8B8pdRcAL8H8H+NykM0EW29LogAU6O0DAHAtBwrJ0NERDFgYUUerpmViy276+F0czo0UU/tP4bTQx6sr7bpjkJxwMjJ0AIA9UqpRqWUG8DzAG4Z+QCl1G6l1FDwyzcAlBiYh2jcTjqcKMhIgcUUvStKp2fzWkNERLFARLB+1Wx09g/j12806Y4Tk3qHPNi6rxErLy7C5aU5uuNQHDDyJ7xiACPnwC3B287nLgB/NzAP0bid7HWdOaQgWk3LtqLN4YRSvLo5EVG0WzArF4ur8vH4ngYMDHt1x4k5v9zfiH6XF+s4FaIwiYpfd4vI1wDMB/Bf57l/tYgcFpHDnZ2dkQ1HCa2t13nm+OpoNS0nFcNeP04P8cKrRESxYP2q2Tg95MFT+4/pjhJTegbdeHL/MXxm7jTMmZ6lOw7FCSPLUCuA0hFflwRv+ycishLAv+I6lmsAABw1SURBVAG4WSk1fK4XUkptVUrNV0rNLyjgRbUoMpRSOOlwRe1JciFnrjXk4L4hIqJYcEVpDlZeXIit+xrRy19kjdkTrzbA6fFh7coq3VEojhhZhg4BqBKRWSKSDOB2AC+PfICIXAngCQSKUIeBWYjGrdfpgdPji9prDIWcudYQ9w0REcWMtdU29Lu8+OX+Rt1RYkJHnwvPvN6EW68oRmVhpu44FEcMK0NKKS+ANQC2ATgK4AWl1BEReVBEbg4+7L8AZAB4UUTeFZGXz/NyRBEXOq46ViZDPFGOiCh2XDI9G5+5bBqe3H8MPYNu3XGi3mN7GuDxKdzPqRCFmdnIF1dKvQLglbNu+8mIz1ca+f5EkxEqF9OifDKUn5ECi0l4rSEiohiztroKf//wJJ54tQE//vTFuuNErTaHE785eAJfvKoEM/PSdcehOBMVBygQRaO24LKz6VE+GUpKEhRl8VpDRESxprIwE7deUYxnXm9CRz9/oXU+j+yqBwDcu4JTIQo/liGi8zjpcMKcJCjITNEdZVTTs1NxkpMhIqKYc9+KKnh8Co/tbtAdJSod7x7Ei4ebcceCUhRH+aUuKDaxDBGdx8leF4qyrDAlie4oo5qWY0UbJ0NERDGnLD8dX7yqBL85eIKngp7D5v+/vTsPk6q+8z3+/nb1TnfTdLM3W7N0I1ECXETFiEgU1EncR8ncjLuOW+bGJRMz3kkmeWai5olxeXRivHNNYoxbyDjBiCGAEJ2JCHgRFLAXEKFZutkaaaDX+t0/+hQp2wYauqpOVZ3P63nq6VPnnDr96V9V16++dc75nSU1hDKMO88b63cUSVMqhkSOYnvjYYYk+TWGIob0zaP+02bCYV14VUQk1UQO/3pyaa3PSZJLbUMT/7l6G9eeNZKBRanRH0vqUTEkchQ79jcfGbY62Q0tzqWtw7G7qdtLdYmISBIrK87ja9OG88rKrWzZc8jvOEnjscXV5GaFuO3cMX5HkTSmYkikG+GwY+f+5iPDVie7yPDf23WtIRGRlHTneWMJZRiPL6nxO0pS2LDjU36/dgc3nl1OaUHyn7srqUvFkEg39hxspbUjnEKHyXnXGtLx5iIiKWlgUS7XTR/Fq6vrqG1o8juO736yqJrC3ExuOWe031EkzakYEunGX64xlCqHyWnPkIhIqvu7GaPJzQrx2OJqv6P4as3WRhatr+fWc0bTNz/L7ziS5lQMiXQjcgHTZL/GUES//CxyMjO0Z0hEJIWVFuRw49nl/H7tDtZv/9TvOL55ZFE1/fKzuOFL5X5HkQBQMSTSjb/sGUqNw+TMjKHFeezQniERkZR2yzmjKcrN5CeLgrl3aMXHe3mrehe3zxxDQU6m33EkAFQMiXRjx/5msjMzKO2T7XeUHhvSV9caEhFJdX3zs7h1xmgWb6hn9ZZ9fsdJKOccP15YxcDCHP72zFF+x5GAUDEk0o3INYbMkv+CqxFD+ubpgn0iImng+rPLKemTHbi9Q2/X7GbF5r3cNWssedkhv+NIQKgYEunG1r2HGN4v3+8YJ2R4SR4NB1o43NrhdxQREemFgpxM7pg5hrdrdrN80x6/4ySEc45H/lhFWXEe15w+3O84EiAqhkS6CIcd1fVNjBtU4HeUE1I5qBDn0JCsIiJp4OtnjmRQUQ4/XliFc87vOHG3aH09a+r28/dfHktOpvYKSeKoGBLpom7fYQ63dVA5qNDvKCdknJe3qv6Az0lERKS3crNC3DVrHKs+2cey6l1+x4mrcNjxk0XVjCrN58opw/yOIwGjYkiki0gxUTE4tYqhUaX5ZIcyqFExJCKSFq6ZOpxh/fJ45I/pvXfotbXb+WjnAe6+oILMkD6aSmLpFSfSRbVXTIwbmFqHyWWGMhgzsEB7hkRE0kR2ZgbfPL+CD7d9yh8+3Ol3nLho6wjz6KJqxg8u5KsTh/odRwJIxZBIF9X1BygrzqMwN/Wuel05qIDqnSqGRETSxeWTyxg7sIBHFlXTEU6/vUO/fa+OzXsOce/sSjIyUmcEV0kfKoZEuqjaeYDKFDtELqJicCHb9zfzaXOb31FERCQGQhnGvRdUUNvQxH+u3uZ3nJhqbuvg8SU1TBpezPmnDPQ7jgSUiiGRKG0dYTbtOphyI8lFRAZ90HlDIiLp48JTB3NqWRGPLq6mtT3sd5yYeeHdLezY38y35lSm1HX9JL2oGBKJ8smeg7R2hFNuJLmIisiIcjs1vLaISLowM+6bXUndvsO8vHKL33Fi4mBLO08trWX6mFLOHtvf7zgSYCqGRKJEioiKFC2Gyorz6JMdOjIIhIiIpIdzKwYwbVQJT7xZmxYX1372vz5mz8FW7ptT6XcUCTgVQyJRquoPkGEwNsVGkovIyDDGDSpUMSQikmbMjH+4sJJdB1r4+Z8/9jtOr+w92Mozb21i9oRBTBnRz+84EnAqhkSiVO88wKjSPuRmpe7VrysGFagYEhFJQ1NHlXD+KQP56bKNNB5q9TvOSfvpsloOtrZrr5AkBRVDIlGqGw6k7CFyERWDCtnd1Mrupha/o4iISIzdN6eSppZ2fvqnjX5HOSnbGw/zy3c+4Yopw1K+v5X0oGJIxNPc1sHm3QepSNGR5CIiw4Jr75CISPoZP7iIyyeV8Yv/3szO/c1+xzlhjy+uAQd3X1DhdxQRQMWQyBEbdzURdp3X6kllfxleWyPKiYiko7svqCDsHI8vqfY7ygmpbWjiN+9t5W/PGklZcZ7fcUQAFUMiR0T2pKTqsNoRAwpzKM7Pokp7hkRE0tLwknz+5xkjeWVVHbUNqfNe//AfPiI/O5M7Zo7xO4rIESqGRDxVO5vIChmj+vfxO0qvmBkVgwqp3pk6HaSIiJyYb8waS15WiIfeqPI7So+s+Hgvi9bXc/vMMZQW5PgdR+QIFUMinur6A4wZUEBWKPX/LSoHFVJVfwDnnN9RREQkDkoLcrh95hgWb6jn3U17/I5zTM45frhgA4OLcrnx7HK/44h8Rup/6hOJkaqdqT+SXETF4EIONLezIwVPrhURkZ658exyBhfl8sMFG5L6y6/XP9jB+1sbuWd2BXnZqXvpCklPKoZEgG2Nh9nWeJiJw/r6HSUmJpZ1/h0rN+/1OYmIiMRLXnaIe2dXsKZuP79fu8PvON1qae/gR3+oYvzgQq6cMszvOCKfo2JIBFhW1QDAzMoBPieJjdPK+lLSJ5tlVbv8jiIiInF0xZRhjB9cyI8WfkRLe4ffcT7n+eVb2LL3EPdfNJ5QhvkdR+RzVAyJAMuqdlFWnMeYAal9jaGIjAxjxrj+vFW9i3A4eQ+dEBGR3gllGP948Sls3XuYZ/9rs99xPmNPUwuPL67mnHH9ObciPb5slPSjYkgCr7U9zJ9rdzOzcgBm6fOt1czKgew52MoH2/b7HUVEROJoRsUAzj9lEE++WUP9p8lzrugji6o52NrBd78yIa36V0kvKoYk8FZt3svB1g5mVg70O0pMzagYgBk6VE5EJAD+6Sun0NbhePiNj/yOAsCH2/bz4ootXHvWSMalyeBEkp5UDEngLaveRXYog+ljSv2OElMlfbKZOKyYpd75UCIikr5GlvbhpnPK+Y/V23jvk32+ZnHO8f3X1tEvP5tvnl/haxaR41ExJIG3rKqB08v70Scn0+8oMXde5QDW1DWy92Cr31FEfGVmF5pZlZnVmtn93SwfaWZLzGytmS0zs2He/Elm9o6ZrfOWXZP49CI9c9d5YxlUlMP3X1vn6/mir63dwcrN+/jWnEr65mX5lkOkJ1QMSaBtbzxMdX0TMyvS6xC5iJmVA3EO3q7RoXISXGYWAp4CLgImAF8zswldVvsx8JxzbiLwA+BBb/4h4Frn3BeAC4HHzKw4MclFTkyfnEzuv2g8a+v28+LKLb5kONDcxg9f38AXhhZx9dThvmQQOREqhiTQIufTpMuQ2l1N1BDbIgDTgFrn3CbnXCvwEnBpl3UmAG9600sjy51z1c65Gm96O9AApOcbhqSFyyaVcdboUh5a8BE7fbjw9sN/+Ij6A838y2WnaihtSQkqhiTQllU1UFacx9iB6TGkdlcaYlsEgDJga9T9Om9etDXAFd705UChmX3mREIzmwZkAxvjlFOk18yMB684jdaOMP/0uw9xLnHv/Ss37+X55Vu4YXo5k0f0S9jvFekNFUMSWK3tYf67djfnptmQ2l1piG2RHrkPONfMVgPnAtuAI1ewNLMhwK+AG5xz4a4PNrNbzWyVma3atUt7YsVfo/r34e4LKli0vp43PtyZkN/Z3NbBt3+7lmH98rhvjgZNkNShYkgCa8EHOzjY2sH5p6Tn+UIRMyoGkBUyXl619fgri6SnbUD0yQvDvHlHOOe2O+eucM5NBh7w5jUCmFkR8DrwgHNueXe/wDn3jHNuqnNu6oABOopO/Hfzl8o5tayI7/5uHfsPtcX99z21tJZNuw7yw8tPIz87/QYkkvSlYkgCqa0jzKOLqzllSFHaDp4QUdInm7mnj+CVlVv5ZM9Bv+OI+GElMM7Mys0sG5gLzI9ewcz6m1mkT/wO8Kw3Pxt4lc7BFeYlMLNIr2SGMnj4yonsO9Qa98PlVm/Zx0+XbeSKKWXMqNCXAZJaVAxJIP1mVR2f7DnEfbMryAjACZ7fmDWWzJDx2OIav6OIJJxzrh24C1gIbABecc6tM7MfmNkl3mozgSozqwYGAf/qzb8amAFcb2bve7dJif0LRE7OF4b25e7zxzF/zXZeWBGf0eUaD7Vy1wurGVSUy/e+8oW4/A6ReNJ+TAmc5rYOnlhSw5QRxcwan957hSIGFuVy3fRRPPPWJm6fOYYKXQ1cAsY5twBY0GXed6Om5wGf2/PjnHseeD7uAUXi5I6ZY1mxeR/ff209XxxWzKllfWO27XDYce8ra2g40My826bTN1/XFJLUoz1DEjjPL/+EnZ82860549N64ISubpsxhoLsTB75Y5XfUUREJEEyMoxHr/4iJfnZ3PnC/+PT5tidP/R/3t7Eko8a+N9/NYEvDtfltyQ1qRiSQGlqaefflm3kS2P7c9aY0uM/II3065PNzeeMZuG6etZsbfQ7joiIJEhpQQ5P/s1k6vYd5p6X19DW8bkBEU/Y2zW7+NHCKv7qtCFce9bIGKQU8YeKIQmMlvYO7nn5ffYebOW+OZV+x/HFTeeUU9Inm7tfed+Xi/GJiIg/po4q4XtfncDiDfV88+X3ae9FQfTn2t3c/MtVjBtYwINXnhaooywk/agYkkA43NrBrc+9xx/X1/O9r05gUkB35xfkZPL01/8H9fubufpn77B17yG/I4mISIJce9YoHrj4FF5fu4N7XllDx0lcjHv5pj3c9MtVjCzN59c3n0FRrs4TktQW12LIzC40syozqzWz+7tZnmNmL3vL3zWzUfHMI8HU1NLO9T9fwVs1u3joitO44exyvyP5alp5Cb++5UwaD7Vy9c/eYdOuJr8jiYhIgtwyYzTfvnA889ds555X3udQa3uPH/un6l3c+IuVlPXL49c3n0lpQU4ck4okRtyKITMLAU8BFwETgK+Z2YQuq90E7HPOjQUeBR6OVx4Jnt1NLTyxpIbzfryMVZ/s47FrJjF32gi/YyWFScOLeenWs2htD3PxE2/zj69+QG3DAb9jiYhIAtw+cwzfmlPJ797fzuxH3+Kt6l3HXL/xUCvf+s0arnt2BcP65fHCLWcwoFCFkKSHeA6tPQ2odc5tAjCzl4BLgfVR61wK/LM3PQ940szMxfPKYJJ2nHM0t4XZc7CFmvom1u/4lA+37WfJRw20toeZWTmAO2aOZVp5id9Rk8qEoUW8esfZPLW0lnnv1fHCu1s4a3QpU0YWM2FIXyoHFzKwKIfCnEwdDy4ikmbuPG8sp48q4f7fruXaZ1dw+eQyLpk0lMnDiynOz6a9I0x1fRPvfryHp5bWsu9QG3fMHMPff3kcuVkhv+OLxEw8i6EyYGvU/TrgjKOt45xrN7P9QCmwOx6BnntnM6+u3haPTUuMRcph591xQNg5wmHoCDvaOsK0doRpaQ+z/3Abre2fPRF0WL88rp46jOunlzN2YEGi46eMEaX5PHzVRP7hwkpeXLGF19bs4Ok/bfrMceShDKNvXhZ5WSGyMzPIChmhjAxCGZBhhgFEfnZOSgqoHFTIQ1dO9DuGiPhoWnkJC/7XOTy1tJan/7TxyGekESX57G5q4VBrBwATh/XluRvPYMLQIj/jisRFSlx01cxuBW4FGDHi5A9zyg5lUJCTEn+ywJG9EQZkmPfB24xQBmRnhsgKGTmZIfrmZVGcn0VxXhajBxQwfkihTug8QaUFOdw1axx3zRpHc1sHtQ1NVNcfYO/BVhoPtdF4uJXmtjCt7WHaOsK0dTicc50FqvOKVjr30klqyMvWN7siArlZIe6dXclt545hTV0jq7c08uG2/QwqymXyiGImD+/H8JI8HSEgaSuelcE2YHjU/WHevO7WqTOzTKAvsKfrhpxzzwDPAEydOvWkP23NnTZC54yIHEduVohTy/rG9CrlIiKS3PrkZDJ9TH+mj+nvdxSRhIrnaHIrgXFmVm5m2cBcYH6XdeYD13nTVwFv6nwhERERERFJhLjtGfLOAboLWAiEgGedc+vM7AfAKufcfOD/Ar8ys1pgL50Fk4iIiIiISNzF9QQa59wCYEGXed+Nmm4G/jqeGURERERERLoT14uuioiIiIiIJCsVQyIiIiIiEkgqhkREREREJJBUDImIiIiISCCpGBIRERERkUBSMSQiIiIiIoGkYkhERERERAJJxZCIiIiIiASSiiEREREREQkkFUMiIiIiIhJI5pzzO8MJMbNdwCe92ER/YHeM4sRLsmdUvt5L9ozK13vJnrG3+UY65wbEKkw6UT+VFJSv95I9o/L1XrJnTEg/lXLFUG+Z2Srn3FS/cxxLsmdUvt5L9ozK13vJnjHZ8wVZKjw3yZ5R+Xov2TMqX+8le8ZE5dNhciIiIiIiEkgqhkREREREJJCCWAw943eAHkj2jMrXe8meUfl6L9kzJnu+IEuF5ybZMypf7yV7RuXrvWTPmJB8gTtnSEREREREBIK5Z0hERERERCQ9iyEz+2szW2dmYTM76igUZnahmVWZWa2Z3R81v9zM3vXmv2xm2THOV2Jmi8ysxvvZr5t1zjOz96NuzWZ2mbfsF2b2cdSySbHM19OM3nodUTnmR81PhjacZGbveK+FtWZ2TdSyuLTh0V5TUctzvPao9dpnVNSy73jzq8xsTizynES+e8xsvddeS8xsZNSybp9rHzJeb2a7orLcHLXsOu81UWNm1/mU79GobNVm1hi1LO5taGbPmlmDmX14lOVmZk94+dea2ZSoZXFvP+mkfioxGb311E999ncmdT/Vw4y+9lXqp3qdL7n6Kedc2t2AU4BKYBkw9SjrhICNwGggG1gDTPCWvQLM9aafBm6Pcb4fAfd70/cDDx9n/RJgL5Dv3f8FcFWc27BHGYGmo8z3vQ2BCmCcNz0U2AEUx6sNj/WailrnDuBpb3ou8LI3PcFbPwco97YT8iHfeVGvs9sj+Y71XPuQ8XrgyW4eWwJs8n7286b7JTpfl/W/ATyb4DacAUwBPjzK8ouBNwADzgTeTVT76faZ50H9VIIyHu3/LhnaEPVTJ5vRt76qh/muR/3UsTImVT+VlnuGnHMbnHNVx1ltGlDrnNvknGsFXgIuNTMDZgHzvPV+CVwW44iXetvt6favAt5wzh2KcY5jOdGMRyRLGzrnqp1zNd70dqABiOdFIrt9TXVZJzr3PODLXntdCrzknGtxzn0M1HrbS2g+59zSqNfZcmBYjDP0OuMxzAEWOef2Ouf2AYuAC33O9zXgxRhnOCbn3Ft0fig9mkuB51yn5UCxmQ0hMe0nHvVTMaF+6sQlez/Vo4w+91Xqp3op2fqptCyGeqgM2Bp1v86bVwo0Oufau8yPpUHOuR3e9E5g0HHWn8vnX6j/6u06fNTMcmKcD3qeMdfMVpnZ8sjhESRhG5rZNDq/IdkYNTvWbXi011S363jts5/O9urJYxORL9pNdH4zE9Hdcx1rPc14pffczTOz4Sf42ETkwztsoxx4M2p2ItrweI72NySi/eTEqJ86NvVTJy7Z+6meZoyW6L5K/VT8JbSfyuztBvxiZouBwd0sesA597tE5+nqWPmi7zjnnJkddUg/rxI+DVgYNfs7dL6xZtM57OC3gR/4lHGkc26bmY0G3jSzD+h84+y1GLfhr4DrnHNhb3ZM2jBdmdnXganAuVGzP/dcO+c2dr+FuHoNeNE512Jmf0fnN5izfMhxPHOBec65jqh5ydKGkgDqp9RP9TCf+qmTlMR9lfqpFJKyxZBz7vxebmIbMDzq/jBv3h46d8dlet+IRObHLJ+Z1ZvZEOfcDu8NsOEYm7oaeNU51xa17cg3TS1m9nPgvhPNF6uMzrlt3s9NZrYMmAz8liRpQzMrAl6n88PH8qhtx6QNuzjaa6q7derMLBPoS+drriePTUQ+zOx8Ojvyc51zLZH5R3muY/0GedyMzrk9UXf/nc7j8iOPndnlscsSnS/KXODO6BkJasPjOdrfkIj2CxT1U4D6KfVT8cnoZ1+lfir+EtpPBfkwuZXAOOscTSabzhfEfOecA5bSefwzwHVArL/Bm+9ttyfb/9yxnN6bauSY58uAbkfjiHdGM+sX2W1vZv2Bs4H1ydKG3vP6Kp3Hnc7rsiwebdjta+oYua8C3vTaaz4w1zpH8SkHxgErYpDphPKZ2WTgZ8AlzrmGqPndPtcxztfTjEOi7l4CbPCmFwKzvaz9gNl89pvqhOTzMo6n8+TOd6LmJaoNj2c+cK11OhPY733oSkT7yYlRP9XLjOqnPifZ+6keZfS5r1I/FX+J7adcnEeM8OMGXE7ncYQtQD2w0Js/FFgQtd7FQDWdFe8DUfNH0/kPXgv8BsiJcb5SYAlQAywGSrz5U4F/j1pvFJ1VcEaXx78JfEDnG+PzQEEc2vC4GYHpXo413s+bkqkNga8DbcD7UbdJ8WzD7l5TdB7WcIk3neu1R63XPqOjHvuA97gq4KI4/W8cL99i738m0l7zj/dc+5DxQWCdl2UpMD7qsTd6bVsL3OBHPu/+PwMPdXlcQtqQzg+lO7zXfh2dx9PfBtzmLTfgKS//B0SNZJaI9tPtSFurn0pAxmP93yVDG6J+6mQz+tpX9SCf+qlj50uqfsq8DYuIiIiIiARKkA+TExERERGRAFMxJCIiIiIigaRiSEREREREAknFkIiIiIiIBJKKIRERERERCSQVQyIiIiIiEkgqhkREREREJJBUDIkkiJmdbmZrzSzXzPqY2TozO9XvXCIiIqB+SoJJF10VSSAz+xc6r+6dB9Q55x70OZKIiMgR6qckaFQMiSSQmWUDK4FmYLpzrsPnSCIiIkeon5Kg0WFyIolVChQAhXR+8yYiIpJM1E9JoGjPkEgCmdl84CWgHBjinLvL50giIiJHqJ+SoMn0O4BIUJjZtUCbc+4FMwsBfzazWc65N/3OJiIion5Kgkh7hkREREREJJB0zpCIiIiIiASSiiEREREREQkkFUMiIiIiIhJIKoZERERERCSQVAyJiIiIiEggqRgSEREREZFAUjEkIiIiIiKBpGJIREREREQC6f8DgUpPngbxPzoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f = lambda x: np.exp(-100*x*x)\n", "a = -1.\n", "b = 1.\n", "L = b-a\n", "N = 99\n", "alpha = 1.0\n", "beta = 1.0\n", "\n", "h=L/(N+1)\n", "# csc_matrix is a sparse matrix data structure that Python claims is more efficient to be used with spsolve.\n", "A = csc_matrix(sp.sparse.spdiags(np.array([np.ones(N)*-2., np.ones(N), np.ones(N)]), [0,-1,1], N,N)) / h**2\n", "\n", "x = np.arange(a,a+L+h,h)\n", "b = f(x[1:-1])\n", "b[0] -= alpha/h**2 # During the lecture, I forgot to divide alpha and beta by h**2\n", "b[-1] -= beta/h**2 # I think I also used + instead of -.\n", "u = np.zeros(N+2)\n", "u[1:N+1] = spsolve(A,b)\n", "u[0] = alpha\n", "u[-1] = beta\n", "fig, ax = plt.subplots(1,2, figsize=(14,6))\n", "ax[0].plot(x,f(x))\n", "ax[0].set_xlabel('x')\n", "ax[0].set_ylabel('f(x)')\n", "ax[1].plot(x,u)\n", "ax[1].set_xlabel('x')\n", "ax[1].set_ylabel('u')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }