jupyter/HW2-AnOddChargeDistribution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerically integrating functions of one variable\n",
    "\n",
    "In this problem, you will learn how to approximately determine the integral of a functions of a single variable using the trapezoidal rule and Simpson's rule. You will compare your results to the simplest built-in integrator from the ``scipy.integrate`` library, ``quad``. In this notebook, we have started by walking through the results of a computing numerically the integral of a function with known anti-derivative, but have left developing the Simpson alogorithm to you."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import all the libraries\n",
    "\n",
    "We will start by importing all the libarries and functions that we need. These can be called upon in any cell, but it's often useful to call them up in the first cell, so you know what's available. Also, to avoid conflicts in namespaces (i.e., functions and methods that appear in multiple libraries), we will use the ``as`` methods of importing libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sympy as sp                 ## Import symbolic python library\n",
    "import numpy as np                 ## Import basic scientific computing library\n",
    "import matplotlib.pyplot as plt    ## Import plotting library\n",
    "from scipy.integrate import quad   ## Import a specific integrator from scipy.intergrate library\n",
    "\n",
    "## Make sure our plots show up inline\n",
    "%matplotlib inline\n",
    "\n",
    "## Add some style to the printing of sympy functions and matplotlib plots\n",
    "sp.init_printing(use_unicode=False, wrap_line=False, no_global=True)\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example\n",
    "\n",
    "We will work with a function with a known anti-derivative,\n",
    "\n",
    "$$ g(x) = 10*e^{-x},$$\n",
    "\n",
    "which are planning to integrate from 0 to 10,\n",
    "\n",
    "$$I(0,10) = \\int_0^{10} 10*e^{-x} dx.$$\n",
    "\n",
    "If this all looks familar, please skip to the problem to solve below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Plot this function\n",
    "\n",
    "First, we will get a sense of what the function looks like over the domain we care about."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110ece748>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAFpCAYAAAB54yVXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VPWh//HPd2aSkBCyDlsICIEAsghGdhdEolWr1tqW\ntlxtrb2L2Nbaert5vbbe/uyTtpfCwy3+aJ/2Yvu7t4tXvVRUXCIKVqwiKcgiEDUoFNoQEgIJkO18\nf39MiGEPk8l8J3Per+fJMzNnzvKBL+qH4/ecY6y1VgAAAIAPBVwHAAAAAFyhDAMAAMC3KMMAAADw\nLcowAAAAfIsyDAAAAN+iDAMAAMC3KMMAAADwLcowAAAAfIsyDAAAAN+iDAMAAMC3KMMAAADwrVC8\nD7h37954H1LhcFg1NTVxPy7ii3H2B8bZHxhnf2Cc/cHFOBcUFHR5Xc4MAwAAwLcowwAAAPAtyjAA\nAAB8izIMAAAA36IMAwAAwLcowwAAAPAtyjAAAAB8izIMAAAA36IMAwAAwLfO+QS6hx9+WBUVFcrO\nztbChQslSQ0NDVq0aJH279+v/v3762tf+5oyMzN7PCwAAAAQS+c8M3zllVfqvvvuO2HZihUrNHHi\nRC1ZskQTJ07UihUreiwgAAAA0FPOWYbHjRt3ylnf9evXa/bs2ZKk2bNna/369T2TLgZs7X41rf+j\n6xgAAABIQOecJnE69fX1ys3NlSTl5OSovr7+jOuWl5ervLxcklRWVqZwOBzNIaPWUL5CB//nEfX/\nr+cVSM+I67ERX6FQKO5/vhB/jLM/MM7+wDj7Q6KPc1RluDNjjIwxZ/y+tLRUpaWlHZ9ramq6e8jz\nYgcWSp6nAxVvyIyZENdjI77C4XDc/3wh/hhnf2Cc/YFx9gcX41xQUNDldaO6m0R2drbq6uokSXV1\ndcrKyopmN/ExvFiSZHftdBwEAAAAiSaqMjxlyhStWbNGkrRmzRpNnTo1pqFiyfTLVnBggWwVZRgA\nAAAnOuc0icWLF2vbtm06fPiw7rzzTs2bN08333yzFi1apNWrV3fcWi2RhYovVNu2t1zHAAAAQII5\nZxm+5557Trv8gQceiHmYnpJSPF5Nf3xRtr5OJjvXdRwAAAAkCF88gS6leFzkDVMlAAAA0Ik/ynDR\naCkQkK2qdB0FAAAACcQXZdik9ZGGXMAdJQAAAHACX5RhSTIjRku7KmU9z3UUAAAAJAjflGGNGC0d\naZSq97lOAgAAgAThmzJsRoyWxMM3AAAA8CHflGENLpTS+kjvUYYBAAAQ4ZsybAJB6YJRsru4owQA\nAAAifFOGJcmMKJZ2vyfb0uI6CgAAABKAz8rwGKm1Vdqzy3UUAAAAJABflWGNKJbERXQAAACI8FcZ\nzg1L2blcRAcAAABJPivDxhhpxBjZ97a7jgIAAIAE4KsyLElm1Fipep/soYOuowAAAMAx/5XhorGR\nN+/tcBsEAAAAzvmuDOuCkVIwJPsuUyUAAAD8zndl2KSmScOKmDcMAAAA/5VhSTIjx0pVlbKtra6j\nAAAAwCH/luGWZml3lesoAAAAcMiXZVjtF9ExVQIAAMDffFmGTV5YygtLXEQHAADga74sw1LkFmvc\nUQIAAMDffFuGNXKsVLtftu6A6yQAAABwxLdl2Iy8MPKGecMAAAC+5dsyrKHDpZRU2XcowwAAAH7l\n2zJsQinS8FGy777tOgoAAAAc8W0ZliIX0emD92Rbml1HAQAAgAP+LsOjxkptrdL777iOAgAAAAd8\nXYbVfhGdrWSqBAAAgB/5ugybftnSoELZyq2uowAAAMABX5dhSTLF46R335b1PNdRAAAAEGe+L8Mq\nHi8daZT2vu86CQAAAOLM92XYFI+TJNnKbY6TAAAAIN58X4aVP0DKDUuUYQAAAN/xfRk2xsgUj5Pd\nuVXWWtdxAAAAEEe+L8OSpOJxUn2ttP+vrpMAAAAgjijDkkzxeEnMGwYAAPAbyrAkDR4qZWRK3G8Y\nAADAVyjDkkwgIBWP48wwAACAz1CG25nicVL1Xtn6OtdRAAAAECeU4XZmVOR+w3qHs8MAAAB+QRk+\n7oKRUmoqUyUAAAB8hDLczoRSpBFjZLmIDgAAwDcow52Y0eOl3VWyRxpcRwEAAEAcUIY7MWMmStZK\nOzk7DAAA4AeU4c6KxkihFNkdW1wnAQAAQBxQhjsxKanSyLGyO95yHQUAAABxQBk+iRk7UdqzS7bx\nsOsoAAAA6GGU4ZOYMRdF5g0zVQIAACDpUYZPNqJYSk2T3bHZdRIAAAD0sFB3Nn7qqae0evVqGWM0\ndOhQ3XXXXUpNTY1VNidMKEUadSFlGAAAwAeiPjNcW1urVatWqaysTAsXLpTneVq3bl0sszljxkyU\n/vK+7OF611EAAADQg7o1TcLzPDU3N6utrU3Nzc3Kzc2NVS6nzJiJkTecHQYAAEhqUU+TyMvL0403\n3qgFCxYoNTVVkyZN0qRJk2KZzZ0LRklp6bI7tshMucx1GgAAAPSQqMtwQ0OD1q9fr6VLlyojI0M/\n+clPtHbtWl1xxRUnrFdeXq7y8nJJUllZmcLhcPcSRyEUCp33cevGT1bbO9uc5EV0ohln9D6Msz8w\nzv7AOPtDoo9z1GV48+bNGjBggLKysiRJ06dP186dO08pw6WlpSotLe34XFNTE+0hoxYOh8/7uF7R\nGNmK17T/3UqZ7OSY/pHsohln9D6Msz8wzv7AOPuDi3EuKCjo8rpRzxkOh8OqrKxUU1OTrLXavHmz\nhgwZEu3uEo4ZHZk3zF0lAAAAklfUZ4aLi4s1Y8YMfetb31IwGNTw4cNPOAPc6w0rktL7StvfkqZd\nce71AQAA0Ot06z7D8+bN07x582KVJaGYYFAaM0H27U2uowAAAKCH8AS6szDjJks1f5Ot3uc6CgAA\nAHoAZfgszIWRW8XZbRsdJwEAAEBPoAyfzcAhUl5Y9m3KMAAAQDKiDJ+FMUbmwsnS9rdkvTbXcQAA\nABBjlOFzGTdZOtIo7XrHdRIAAADEGGX4HJg3DAAAkLwow+dg+mVLw4qYNwwAAJCEKMNdYC6cLL27\nQ/bYUddRAAAAEEOU4S4w4yZLba1S5VbXUQAAABBDlOGuKB4npaQybxgAACDJUIa7wKSkSsXjKMMA\nAABJhjLcRWbcZGnvB7IHD7iOAgAAgBihDHeRGXexJMlu5ewwAABAsqAMd1XhcCk7T9qywXUSAAAA\nxAhluIuMMTITLpbd9mfZNh7NDAAAkAwow+fBTLgk8mjmqh2uowAAACAGKMPnY9xkKRCQ3cxUCQAA\ngGRAGT4PJiNTKhory7xhAACApEAZPk9mQon0wXuy9XWuowAAAKCbKMPnyUy8RJJkt1Q4TgIAAIDu\nogyfr6FFUnYut1gDAABIApTh8xS5xVoJt1gDAABIApThKHCLNQAAgORAGY4Gt1gDAABICpThKHCL\nNQAAgORAGY6SmXhJ5BZrdQdcRwEAAECUKMNRMpOmS5LsW+sdJwEAAEC0KMPRKhgqhQfKbnrDdRIA\nAABEiTIcJWOMzKRp0tubZJuOuY4DAACAKFCGu8FMmia1tkjbNrqOAgAAgChQhrujeLyU3ld20+uu\nkwAAACAKlOFuMKFQ5Gl0b70p6/E0OgAAgN6GMtxdk6ZJh+ulqkrXSQAAAHCeKMPdZCZcIgWDTJUA\nAADohSjD3WT6Zkqjxslu5BZrAAAAvQ1lOAbM5GnSvt2y1ftcRwEAAMB5oAzHgLlomiTxAA4AAIBe\nhjIcA2bAYKlgmOxG5g0DAAD0JpThGDElM6XKbbKHDrqOAgAAgC6iDMeIKZklWY+zwwAAAL0IZThW\nCodL/QfJVqxznQQAAABdRBmOEWNMZKrE9rdkjzS4jgMAAIAuoAzHkLl4ptTWJrtpvesoAAAA6ALK\ncCyNGC3l5MtWvOY6CQAAALqAMhxDJhCQuXiGtLVC9thR13EAAABwDpThGDOXzJJamqWtFa6jAAAA\n4Bwow7E2apyUmSW7gbtKAAAAJDrKcIyZYFDm4hmyb70p29LsOg4AAADOgjLcA8zFM6Wmo9K2ja6j\nAAAA4Cwowz3hwoukvv1k17/iOgkAAADOItSdjRsbG7Vs2TLt3r1bxhgtWLBAo0ePjlW2XsuEUmRK\nZsq+8Ypsc5NMaprrSAAAADiNbpXh5cuXa/Lkybr33nvV2tqqpqamWOXq9czUy2VfeV7avEG6ZJbr\nOAAAADiNqKdJHDlyRG+//bauuuoqSVIoFFLfvn1jFqzXGz1B6pctb/1a10kAAABwBlGfGa6urlZW\nVpYefvhhvf/++yoqKtLtt9+uPn36xDJfr2WCQZkpl8r+sVz22BGZPhmuIwEAAOAkUZfhtrY2VVVV\n6Y477lBxcbGWL1+uFStW6DOf+cwJ65WXl6u8vFySVFZWpnA43L3EUQiFQk6O21x6o+peekaZ721X\n+hXXxP34fuNqnBFfjLM/MM7+wDj7Q6KPc9RlOD8/X/n5+SouLpYkzZgxQytWrDhlvdLSUpWWlnZ8\nrqmpifaQUQuHw06Oa8ODpZx8HVr9jBrHlcT9+H7japwRX4yzPzDO/sA4+4OLcS4oKOjyulHPGc7J\nyVF+fr727t0rSdq8ebMKCwuj3V1SMoGAzNTLpC0Vso0NruMAAADgJN26m8Qdd9yhJUuWqLW1VQMG\nDNBdd90Vq1xJw0y9XPaFP8hu/JPMpaXn3gAAAABx060yPHz4cJWVlcUqS3IaXiyFB8q+8YpEGQYA\nAEgoPIGuhxljZKZeLm3fJHuoznUcAAAAdEIZjgMzc47kebJvcM9hAACAREIZjgMzeKh0wSjZ1152\nHQUAAACdUIbjxMy4UvrgXdm9H7iOAgAAgHaU4Tgx066QAgHZP73kOgoAAADaUYbjxGTlSONLZF9f\nI+t5ruMAAABAlOG4MjOulGprpJ1bXEcBAACAKMNxZSZPl/qkM1UCAAAgQVCG48ikpslcMkt2wzrZ\npibXcQAAAHyPMhxnZsYc6dhR2Y1/ch0FAADA9yjD8TZ6gpTXX/a11a6TAAAA+B5lOM5MICAza660\nbaPsgf2u4wAAAPgaZdgBc+lcSZJd96LjJAAAAP5GGXbAhAdKYy+SfbWcew4DAAA4RBl2xFxaKh2o\nlra/5ToKAACAb1GGHTElM6WMTNlXy11HAQAA8C3KsCMmJVVm+mzZitdkGxtcxwEAAPAlyrBD5rKr\npdYW2ddfdh0FAADAlyjDDplhRdKwItk/vuA6CgAAgC9Rhh0zl10t7a6Sff8d11EAAAB8hzLsmJk+\nW0pNk13zrOsoAAAAvkMZdsxkZMpMu0L29TWyRxpdxwEAAPAVynACMLOvlZqbZP/0kusoAAAAvkIZ\nTgBmeLF0wSjZl1fJWus6DgAAgG9QhhOEufI6ad9uqXKr6ygAAAC+QRlOEGbqFVJGX9mXV7mOAgAA\n4BuU4QRh0tJkZl4VeSLdoTrXcQAAAHyBMpxAzOzrpLZW2VdfdB0FAADAFyjDCcQMLpTGTJRd86ys\n1+Y6DgAAQNKjDCeYwFUflQ5USxvfcB0FAAAg6VGGE83k6VL+AHkvPuk6CQAAQNKjDCcYEwjKXPVR\naedW2Q/ecx0HAAAgqVGGE5C57GoprY/siytdRwEAAEhqlOEEZDIyI7dZe2MNt1kDAADoQZThBGXm\n3iC1tsquec51FAAAgKRFGU5QZlChNOES2TWrZFtaXMcBAABISpThBBYovUmqr5N984+uowAAACQl\nynAiGzdZGjxU9vn/lbXWdRoAAICkQxlOYMYYmY/cIu3ZJW39s+s4AAAASYcynODM9CuknHx5zz3h\nOgoAAEDSoQwnOBNKkbn6Jmn7W7K7Kl3HAQAASCqU4V7AXP4RKb2v7LOcHQYAAIglynAvYNIzZK68\nTrbiNdnqva7jAAAAJA3KcC9h5t4oBQOyL/zBdRQAAICkQRnuJUx2buQRza++yCOaAQAAYoQy3IuY\nj9wSeUTz85wdBgAAiAXKcC9iBhbITL1c9uVnZA8fch0HAACg16MM9zLmo5+Smptky590HQUAAKDX\nowz3MqZgmEzJLNnVK2UbG1zHAQAA6NUow72QuWGedOyo7IsrXUcBAADo1bpdhj3P0ze/+U2VlZXF\nIg+6wBSOkCbPkH3xSdmjR1zHAQAA6LW6XYafeeYZDRkyJBZZcB4CN3xaOtIou/op11EAAAB6rW6V\n4QMHDqiiokJz586NVR50kblgpHTRVNkX/iB7pNF1HAAAgF6pW2X4kUce0a233ipjTKzy4DwEbpov\nNR6WLee+wwAAANEIRbvhhg0blJ2draKiIm3duvWM65WXl6u8vFySVFZWpnA4HO0hoxYKhZwct8eF\nwzo4c46aX3hSeZ/8nAJZOa4TOZW044wTMM7+wDj7A+PsD4k+zsZaa6PZ8De/+Y3Wrl2rYDCo5uZm\nHT16VNOmTdPdd9991u327t0bVdDuCIfDqqmpiftx48Hu2y3vu1+RufpjCnzqC67jOJXM44wPMc7+\nwDj7A+PsDy7GuaCgoMvrRn1meP78+Zo/f74kaevWrVq5cuU5izBizwweKjNjtuxLT8tefZNMTr7r\nSAAAAL0G9xlOAubGz0pem+zTj7qOAgAA0KtEfWa4s/Hjx2v8+PGx2BWiYPoPkrn8GtlXnpe95uMy\n/Qe5jgQAANArcGY4SZiPzpMCQdkV/+06CgAAQK9BGU4SJidf5uqPyb6xRnZXpes4AAAAvQJlOImY\naz8h9cuW9z/LFeVNQgAAAHyFMpxETHqGzE2flXZukTa97joOAABAwqMMJxlz+UekQYXyHvuVbGur\n6zgAAAAJjTKcZEwwqMAnvyD97S+yrzznOg4AAEBCowwno4umSGMmyj75W9kjja7TAAAAJCzKcBIy\nxijwqTukhkOyqx5zHQcAACBhUYaTlLlgpMzMq2TL/yD7t72u4wAAACQkynASM5/4vBRKkfe7n3Or\nNQAAgNOgDCcxk50rc9N8aUsFt1oDAAA4DcpwkjNzPioVDJP3u1/INje5jgMAAJBQKMNJzoRCCsz/\nJ+lAteyzj7uOAwAAkFAowz5gxkyUmXq57KrHZff/1XUcAACAhEEZ9gnzyS9IwaC83//CdRQAAICE\nQRn2CZMXlrnh09KmN2QrXnMdBwAAICFQhn3ElH5MKhwh7zc/kz3S4DoOAACAc5RhHzGhkAKf/7J0\n6KDs4792HQcAAMA5yrDPmOHFMqU3yq59VnbnFtdxAAAAnKIM+5D52N9J+QPk/b+lsi3NruMAAAA4\nQxn2IZPWR4HbviT99S+yTz/qOg4AAIAzlGGfMuMvlpk5R/bZx2X3VLmOAwAA4ARl2MfMvC9KffvJ\n++Vi2dYW13EAAADijjLsYyYzS4HPfVnaUyW78neu4wAAAMQdZdjnzKRpMpeWRh7V/O5213EAAADi\nijIMmU//vZQXlvefi2WbjrmOAwAAEDeUYcikZyhw+91S9V7Zx3/lOg4AAEDcUIYhSTJjL5IpvUn2\npadlt/3ZdRwAAIC4oAyjg/n4bdLgoZHpEocOuo4DAADQ4yjD6GBS0xT4x3+WGhvkLV8s63muIwEA\nAPQoyjBOYApHRC6o21Ih+8IfXMcBAADoUZRhnMLMvlYqmSX7v7+WrdrpOg4AAECPoQzjFMYYBT7/\nZSknX97Pfyx7pNF1JAAAgB5BGcZpmYxMBf7hn6Xa/bK//qmsta4jAQAAxBxlGGdkRo6VueVzshte\nZf4wAABISpRhnJW55uOR+cOPPyK7/S3XcQAAAGKKMoyzMsYo8IW7pYFDIvOHa/e7jgQAABAzlGGc\nk+mTocBd35FamuX93zLZlmbXkQAAAGKCMowuMYMKFfjCPdKuStnf/tx1HAAAgJigDKPLTMlMmes/\nJfvK8/JefMp1HAAAgG6jDOO8mI/NlyZPl/39L2Q3b3AdBwAAoFsowzgvJhBU4O/vlYYOl/fzH8n+\n5X3XkQAAAKJGGcZ5M2l9FPjS/VKfdHn/8X3ZQ3WuIwEAAESFMoyomLywAl++Xzp8UN7SH8g2N7mO\nBAAAcN4ow4iauWCUAl/8ulS1U94vFsp6ba4jAQAAnBfKMLrFlMyS+fTfS3/+k+x/L5O11nUkAACA\nLgu5DoDeLzD3Rnn1dbKrHpOyc2Vumu86EgAAQJdQhhET5uO3SYcOyq78nbysHAWuvN51JAAAgHOi\nDCMmjDHSbV+SPVwv+5ufyWZmyUy5zHUsAACAs2LOMGLGBIMK/OM3pZFjIxfUbXzddSQAAICzogwj\npkxamgJfeUAaWiTvZz+U3cJT6gAAQOKKeppETU2Nli5dqoMHD8oYo9LSUl1/PfNEIZmMvgrc86C8\nn9wvb+kPFLj7AZkLJ7mOBQAAcIqozwwHg0HddtttWrRokR566CE999xz2rNnTyyzoRczfTMVuOff\npIEF8n76fdkdW1xHAgAAOEXUZTg3N1dFRUWSpPT0dA0ZMkS1tbUxC4bez/TLUuDr35fyB8pb8qDs\n25tcRwIAADhBTOYMV1dXq6qqSqNGjYrF7pBETFaOAvf+H6n/IHlL/k1285uuIwEAAHQwtpuPDDt2\n7Ji++93v6pZbbtH06dNP+b68vFzl5eWSpLKyMjU3N3fncFEJhUJqbW2N+3HxIe9QveoevEetH7yr\n7K8/qD4z58T8GIyzPzDO/sA4+wPj7A8uxjk1NbXL63arDLe2tuqHP/yhJk2apBtuuKFL2+zduzfa\nw0UtHA6rpqYm7sfFieyRRnlLHpSqdsp84R4FZlwZ0/0zzv7AOPsD4+wPjLM/uBjngoKCLq8b9TQJ\na62WLVumIUOGdLkIw9+O32VCoyfI/ucieS8+5ToSAADwuajL8I4dO7R27Vpt2bJF3/jGN/SNb3xD\nFRUVscyGJGT6pCvwlX+VJk2X/d3P5T3+K3Vzpg4AAEDUor7P8NixY/Xoo4/GMgt8wqSmKbDgW7K/\n/bnss49LB2ulz39ZJpTiOhoAAPCZqMsw0B0mEJTm3ynl5Muu+C/ZQwcVWPAtmT4ZrqMBAAAf4XHM\ncMYYo8BH58ncfre0fZO8sm/J1vzNdSwAAOAjlGE4F7i0VIGvfk+qq5H30L2yO7e6jgQAAHyCMoyE\nYMZNVuA7/y5l9pP3k3+V98rzriMBAAAfoAwjYZhBQxT4zo+lMRNlf/1Teb//hWxbm+tYAAAgiVGG\nkVBMRqYCdz8gM/dG2fIn5S16QLa+znUsAACQpCjDSDgmGFTgM/8g84V7pKod8r5/j+yOLa5jAQCA\nJEQZRsIKzLoqMo+4T4a8hffLW/W4rOe5jgUAAJIIZRgJzRQOV+BfFsqUzJR94lfyHv6BbMMh17EA\nAECSoAwj4Zn0DJl/+qbMZ/5B2lIh73t3y277s+tYAAAgCVCG0SsYYxSYe6MC9/1Yyugrb9F3I3eb\naGl2HQ0AAPRilGH0KmbYSAX+5Scyc66P3G3ioXtl9+xyHQsAAPRSlGH0OiYtTYH5dypw9wPSoYPy\nHvq6vKcflW1tdR0NAAD0MpRh9Fpm4hQFvvcfMpOmy674L9V+84uyH7zrOhYAAOhFKMPo1UxWjgJ3\nfkuBBd+Rd7BW3kP3ynvi18wlBgAAXUIZRlIwJTOVv+S/ZWbOkV31mLwHvyq7baPrWAAAIMFRhpE0\nAplZCtz+VQXueVDy2uQtekDez34kW1vjOhoAAEhQlGEkHTP+YgUe/KnMx/5OdtMb8h64S95zT8i2\ntriOBgAAEgxlGEnJpKQqcMOnFXjwp9LYi2Qfe0Teg3fLbnxd1lrX8QAAQIKgDCOpmf6DFPzy/Qp8\n5V8lGXlLH5L37/fJVlW6jgYAABJAyHUAIB7MRVMVGF8i+8rzsk/+Rt4P7pWZNlvm47fKhAe6jgcA\nAByhDMM3TDAoc+V1sjNmyz77hOwLK2QrXpW54lqZ6z4hk5PvOiIAAIgzyjB8x/TJkLn5VtkrrpV9\n6neya1bJvvK8zOxrZa79hEx2ruuIAAAgTijD8C2TF5b53Jdlr/uk7NO/l139lOzaZ2VmXydz7S0y\nWZRiAACSHWUYvmf6D5K5/auy139K9qnfy5avlH15lcysq2SuuVlmQIHriAAAoIdQhoF2ZkCBzB1f\nk71+nuzz/yv7arns2udlSmZGzhQPL3YdEQAAxBhlGDiJGTQkMn3ipvmyq1fKvvys7IZXpTETFSi9\nUbpoqkwg6DomAACIAcowcAYmJ0/mls/LXvcp2Veeky1fKW/pD6T8AZGL7S67RqZfluuYAACgGyjD\nwDmY9AyZaz4uO/cmadPr8l56RvaJX8s++VuZqZfJzLlBZgRTKAAA6I0ow0AXmWBQKpmlYMks2b0f\nyL70jOxrL8m+9pJUOELm0rky06/kbDEAAL0IZRiIgikYJvN3d8re8jnZP70cudju97+QfewRadI0\nBS6dK40viRRoAACQsCjDQDeY9AyZOddLc66X3bNL9tUXZV9/WV7FOik7T2bq5TLTLpeGF8sY4zou\nAAA4CWUYiBFTOFzm01+U/cTnpM0b5K17Ufblp2XL/yCFB0aK8dTLpcLhFGMAABIEZRiIMRNKkS6e\noeDFM2SPNMj++XXZ9Wtln3tCdtVj0qDCyIV3k2dIQ0dQjAEAcIgyDPQgk5Epc+lc6dK5sofrZSte\nk13/SuRJdyt/J+X1l5k0TWbydGn0+EiRBgAAcUMZBuLE9MuWmX2tNPta2UN1sm+9KbvxddlXX5B9\n6WkpPUNmwiWRh3qMmySTles6MgAASY8yDDhgsnJlLrtauuxq2aYm6e2NspvekN30hrT+FVkpMoVi\n3MUy4y+WRl0ok5LqOjYAAEmHMgw4ZtLSpMnTZSZPl/U8afd7slv/HPkpf1L2uSek1FRp9ASZCyfL\njJkoDR3OI6EBAIgByjCQQEwgIF0wSuaCUdL1n5I9dlTauUV228ZIOf6f/4ycNU7PkEZeKDN6vEzx\neGn4KOYbAwAQBcowkMBMn/TIHOKLpkqSbG2NbOVWqXKr7M6tsk9siJTj1FRpxBiZkRdGHg1dNJo5\nxwAAdAFlGOhFTF5YZvpsafpsSZI9XC9VbpM9Xo6ffSwy1UKS8gfIDG8vxiPGSMNGRqZkAACADpRh\noBcz/bIF9+ggAAAL20lEQVSlkpkyJTMlKXIx3gfvylbtkKoqZat2ShtejZw9DgSkwUNlho6QCkdE\nXocWyfTLcvprAADAJcowkERMWppUPE6meFzHMnuoTnpvp2xVpezu92S3b5b+9HKkIEtSTn7kzhVD\niyJPxysYKg0okElhDjIAIPlRhoEkZ7JyO+5WcZw9fEjaUyW7+z1pd5Xs7irZrRWS5314Fjk8SCoY\nKjO4UBo0NFKSBxVG5jEDAJAkKMOAD5l+WdKFk2QunNSxzLY0S/v2yO7bLe3bLbtvT+R185tSW9uH\nZ5Jzw9KAwTIDBkv9B8sMGCT1HywNGCTTJ8PJrwcAgGhRhgFIUuShHsOKZIYVnbDctrZK+/8q7fsg\nUpD/ukd2/19lN74uHa7/sCRLUr/sSFHuP1jqPzDyuOm8/lJefykvLJPKBXwAgMRCGQZwViYUkgYX\nSoMLZU76zh49Iu3fJ+3/q2z1h692x2bp9Zcla08ty7nhSEnOby/Ief0jy7JzpexcnrQHAIgryjCA\nqJn0DGnYyMht2076zra2SAdrpdr9srX7pQP7pdqayPvqvbLbN0nHjp5YliUpo6+UnRcpxtm57e9z\npOw8mexctQ4vkm31pIzMyENKAADoBsowgB5hQilSeKAUHnhKUT7OHmmUaqululrZ+lqpvk6qr5Wt\nPxh5fXd7ZFlLc2R9SQc6DhCQMvtJmVlSvywpM0smM/ukz1mRs9GZ/aSMTKlPuow5UxoAgB9RhgE4\nYzL6Shnt9z0+wzrWWuloY6QUH6xVpteqw/v+Ih0+JDUckm04JDXURy7+a9gmNRyWbOTBI6ecdTaB\nyJnnjL6RcpzRV0rv254j84TvTPt3Su8r9ekj9UmX0vrw2GsASDLdKsMbN27U8uXL5Xme5s6dq5tv\nvjlWuQBAkiJncjMyIz+Dhyo9HFZjTc0Z17eeFynPh+ulhvbCfPhQZFljo3S0QTrSGDkrfaRBOlgr\ne7RROtIoNTd9uJ8zHSAUai/G6R0FWX0i783x9yd/l5omk5oqpaZJKWmR147PqR3vTSAY2988AMA5\nRV2GPc/TL3/5S91///3Kz8/Xd77zHU2ZMkWFhYWxzAcA58UEAlLffpGf48u6uK1taYmU5iON7eW5\nQTp2RPbYUenYUanpWOT12FGp6Whk+fFlB2tlO3/f1nrivrsSIBg6TVE+8bMJpUQKeSgl8pOSEtnu\n+PuTvw+lRC6CPL48JUUKpkgpndcJRfYRCErBTj8mwLQSAEkv6jL8zjvvaNCgQRo4cKAkadasWVq/\nfj1lGECvZVJSpJQcKSvnxOVR7Mu2tnxYjFuapebmyJnnliapuVm2ual9eeRzZHlTp/Xa12lujqzX\n2CA1N0VuddfaInW8tv+0tZ05SxT5OwRDUjDQ/hpsL8yhyINZji/rWB48w7KgzPFtAsH21/YfEzjp\ns5ECATX0zZR3rOmEZafdxnTe7sR9m87bmaAUMCdtaySZyOu5fs62nkxk36dbp/N3XVnHBCJ/4M6w\nb/5yAsRe1GW4trZW+fn5HZ/z8/NVWVkZk1AA0NuZUIqUmRK5oO9038f4eNbzImejWzoV5M6FuaX9\nc1uL1BJZbjuv57VFtm9r34/XFinYp13W/uO1yba1nrJMba0n7bdNtvO2nheZ19351fMkazveNx5f\n3p3fkxj93iacjkJsPnw53bLjb8xplnXs56TvOu/ntOsf/3jSss7ZziPL/mBQXpt34n7MafZ9yj67\nkOV0zrROtH/JOOt2ZzpWjPd3zu2i2CaG+zMXjJLu/pfz318c9fgFdOXl5SovL5cklZWVKRwO9/Qh\nTxEKhZwcF/HFOPsD4+wPoVBILS0tnQpypETb4wW6LbLMdpToTt97bZJnP/z++DJrO33fvp21kYs0\nrZVkI9tZGzmGlaz9cL0T1pFtf3x5521OXcd63knbtL/a49urvfSfuF9J7dkiOXQ8h3e84tuOdSLH\n7FT+bafvOtY76TvZk9br5PjvyQnH6bT92b7TSfvsnO2U41kFTEDe8bE4+dd1crbj+zrh13C6X/OZ\nnOG7s24Sxf7Ott1Z9nfWv7ydLccZvzvbNlEe64y/h2feJpSTk/D/3o66DOfl5enAgY6bHOnAgQPK\ny8s7Zb3S0lKVlpZ2fK45y4UvPSUcDjs5LuKLcfYHxtkfwuHwCf+NOTsjmZAUVOQHvQb/PCe/Nkmt\nra1xH+eCgoIurxv1HetHjhypffv2qbq6Wq2trVq3bp2mTJkS7e4AAACAuIv6zHAwGNQdd9yhhx56\nSJ7nac6cORo6dGgsswEAAAA9qltzhktKSlRSUhKrLAAAAEBcRT1NAgAAAOjtKMMAAADwLcowAAAA\nfIsyDAAAAN+iDAMAAMC3KMMAAADwLcowAAAAfIsyDAAAAN+iDAMAAMC3KMMAAADwLWOtta5DAAAA\nAC744szwt7/9bdcREAeMsz8wzv7AOPsD4+wPiT7OvijDAAAAwOlQhgEAAOBbwe9973vfcx0iHoqK\nilxHQBwwzv7AOPsD4+wPjLM/JPI4cwEdAAAAfItpEgAAAPCtkOsAPW3jxo1avny5PM/T3LlzdfPN\nN7uOhBirqanR0qVLdfDgQRljVFpaquuvv951LPQAz/P07W9/W3l5eQl/dTKi09jYqGXLlmn37t0y\nxmjBggUaPXq061iIsaeeekqrV6+WMUZDhw7VXXfdpdTUVNex0E0PP/ywKioqlJ2drYULF0qSGhoa\ntGjRIu3fv1/9+/fX1772NWVmZjpOeqKkPjPseZ5++ctf6r777tOiRYv06quvas+ePa5jIcaCwaBu\nu+02LVq0SA899JCee+45xjlJPfPMMxoyZIjrGOhBy5cv1+TJk7V48WL9+Mc/ZryTUG1trVatWqWy\nsjItXLhQnudp3bp1rmMhBq688krdd999JyxbsWKFJk6cqCVLlmjixIlasWKFo3RnltRl+J133tGg\nQYM0cOBAhUIhzZo1S+vXr3cdCzGWm5vbMTE/PT1dQ4YMUW1treNUiLUDBw6ooqJCc+fOdR0FPeTI\nkSN6++23ddVVV0mSQqGQ+vbt6zgVeoLneWpublZbW5uam5uVm5vrOhJiYNy4caec9V2/fr1mz54t\nSZo9e3ZC9rCkniZRW1ur/Pz8js/5+fmqrKx0mAg9rbq6WlVVVRo1apTrKIixRx55RLfeequOHj3q\nOgp6SHV1tbKysvTwww/r/fffV1FRkW6//Xb16dPHdTTEUF5enm688UYtWLBAqampmjRpkiZNmuQ6\nFnpIfX19x192cnJyVF9f7zjRqZL6zDD85dixY1q4cKFuv/12ZWRkuI6DGNqwYYOys7MT+tY86L62\ntjZVVVXpmmuu0Y9+9COlpaUl5P9SRfc0NDRo/fr1Wrp0qX72s5/p2LFjWrt2retYiANjjIwxrmOc\nIqnLcF5eng4cONDx+cCBA8rLy3OYCD2ltbVVCxcu1OWXX67p06e7joMY27Fjh95880196Utf0uLF\ni7VlyxYtWbLEdSzEWH5+vvLz81VcXCxJmjFjhqqqqhynQqxt3rxZAwYMUFZWlkKhkKZPn66dO3e6\njoUekp2drbq6OklSXV2dsrKyHCc6VVKX4ZEjR2rfvn2qrq5Wa2ur1q1bpylTpriOhRiz1mrZsmUa\nMmSIbrjhBtdx0APmz5+vZcuWaenSpbrnnns0YcIE3X333a5jIcZycnKUn5+vvXv3SoqUpsLCQsep\nEGvhcFiVlZVqamqStVabN2/mQskkNmXKFK1Zs0aStGbNGk2dOtVxolMl/UM3Kioq9Ktf/Uqe52nO\nnDm65ZZbXEdCjG3fvl0PPPCAhg0b1vG/Xz772c+qpKTEcTL0hK1bt2rlypXcWi1J7dq1S8uWLVNr\na6sGDBigu+66K+Fuw4Tue/TRR7Vu3ToFg0ENHz5cd955p1JSUlzHQjctXrxY27Zt0+HDh5Wdna15\n8+Zp6tSpWrRokWpqahL21mpJX4YBAACAM0nqaRIAAADA2VCGAQAA4FuUYQAAAPgWZRgAAAC+RRkG\nAACAb1GGAQAA4FuUYQAAAPgWZRgAAAC+9f8BJbrHnukK//QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110cf2588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.arange(0,10,0.01)\n",
    "g = 10*np.exp(-t)\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.plot(t,g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Does it have a known anti-derivative?\n",
    "\n",
    "Next, we can use the ``sympy`` library to check if this function has a known anti-derivative, that is, can it be integated analytically? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEMAAAARBAMAAAB0ogy8AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMlTvq5l2Zoki\nu0Rn3bgMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA6ElEQVQoFWNgwAt2tkx3wKuAgVOgiPkBqhLG\n7tUomlgYtJAVsCYyMDAFMNpAxGLvAoEAA8MfJCXi7R8ZGJYxMGxCEmNgYPnAgsTnBCq5zsAQD9QK\nB+e5DbjhHAYGkJLvDAzvA5DExOft3IjEBSph/AFUMoFBemYPkjgSE6iE5zcDg/8GPhuG00jiSEyQ\nEqAp/huOBTC2IIkzMAgZg4AK2C0Qi8pWzkJ2MpJqqHPjA4BuxgFASu4yMKwXADoIBwApAQZdE0M2\n0FLsakBKmAMYUxgmMTBi9zRH2acqBsadwLjn7FyF3RDcogBgYzyStFtioQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- 10 e^{- x}$$"
      ],
      "text/plain": [
       "     -x\n",
       "-10*e  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=sp.Symbol('x')\n",
    "sp.integrate(10*sp.exp(-x),x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the function has a known anti-derivative, so numerically finding the integral from 0 to 10 is not necessary. We can simply use the antiderivative to evaluate it at the end points and compute the integral,\n",
    "\n",
    "$$I = -10*e^{-x}$$\n",
    "\n",
    "So that,\n",
    "\n",
    "$$I(0,10) = -10({e^{-10}-1}) \\approx 9.9995460007$$\n",
    "\n",
    "But let's see what happens when we numerically integrate it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Define the function\n",
    "\n",
    "First, we define the function, so we can call it as we use our trapezoidal method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def g(t):\n",
    "    return 10*np.exp(-t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Trapezoidal Approach\n",
    "Now we apply the trapezodial algorithm to this integral to determine the approximate value. Here, we are using 10 equally sized steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.819275852\n"
     ]
    }
   ],
   "source": [
    "N = 10\n",
    "a = 0\n",
    "b = 10\n",
    "h = (b-a)/N\n",
    "s = 0       # this just ensures we don't add more to it if we run it again!\n",
    "\n",
    "s = 0.5*g(a) + 0.5*g(b)\n",
    "for k in range(1,N):\n",
    "    s += g(a+k*h)\n",
    "    \n",
    "I = s*h\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our approximate answer isn't great. It's off by about 8%. We can increase the number of steps to get better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.0078775672\n"
     ]
    }
   ],
   "source": [
    "N = 100\n",
    "a = 0\n",
    "b = 10\n",
    "h = (b-a)/N\n",
    "s = 0       # this just ensures we don't add more to it if we run it again!\n",
    "\n",
    "s = 0.5*g(a) + 0.5*g(b)\n",
    "for k in range(1,N):\n",
    "    s += g(a+k*h)\n",
    "    \n",
    "I = s*h\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With 10 times more steps, we are now off by only about 0.08%! But we can do even better with more steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.99962933011\n"
     ]
    }
   ],
   "source": [
    "N = 1000\n",
    "a = 0\n",
    "b = 10\n",
    "h = (b-a)/N\n",
    "s = 0       # this just ensures we don't add more to it if we run it again!\n",
    "\n",
    "s = 0.5*g(a) + 0.5*g(b)\n",
    "for k in range(1,N):\n",
    "    s += g(a+k*h)\n",
    "    \n",
    "I = s*h\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Can we do even better?\n",
    "\n",
    "Now, we are only off by 0.0007%! We can keep doing this for a while, but at some point we are running too many iterations for the gain in precision. We will also run into the issue of machine precision, which is error inherent in the storing on numbers in a computer. \n",
    "\n",
    "There are better methods for numerical integration, which are inherently faster (fewer iterations for same or better precision) and/or adaptive to the problem (get you the best answer for your required precision). One of these is Gaussian quadrature, which the ``quad`` function performs. The first result is the integral result. What's the second one?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9.999546000702377, 2.83275284251709e-13)\n"
     ]
    }
   ],
   "source": [
    "I = quad(g, 0, 10)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problems to Solve\n",
    "\n",
    "Up till now, most of your experience with integration has likely been integrating functions that have anti-derivatives. While this is kind of integration is problematic when you have data that must be integrated (e.g., using measures of position to determine potential energy), we will focus on functions for which there are no anti-derivatives.\n",
    "\n",
    "Consider a line of charge that lives on the x-axis. It exists from $x=0$ to $x=3$, and distribution of that charge is given by,\n",
    "\n",
    "$$\\lambda(x) = 0.5e^{-x^2}.$$\n",
    "\n",
    "In this problem, you will work through this Jupyter notebook to determine the total charge one this line by performing the integral,\n",
    "\n",
    "$$Q = \\int_0^3 0.5e^{-x^2} dx.$$\n",
    "\n",
    "As you work through this problem, you will work through the following activities:\n",
    "\n",
    "1. Plot the distribution of charge. Does it look familiar to you in some way? Maybe you've seen it in other classes?\n",
    "2. Use ``sympy`` to compute the integral of the function (i.e., the total charge). What does it return? Discuss the function. Is it analytic?\n",
    "3. Apply the trapezoidal rule to compute the integral of this function. Use 10 equal width steps. To receive full credit, you must write the code to do this and not use built-in integration functions of ``scipy``.\n",
    "4. Apply Simpson's rule to compute the integral of this function. Again, use 10 equal with steps. Compare the result to part 3. To receive full credit, you must write the code to do this and not use built-in integration functions of ``scipy``.\n",
    "5. Compute both integrals with increasingly more steps (i.e., from 10 to 100) and plot the results. That is compare the plots of what your integration schemes produce on a figure where the number of steps in on the horizontal axes and the value of the integral is on the vertical.\n",
    "6. Look up ``scipy.integrate`` built-in quadrature function, ``quad.`` Use it to compute the same integral and compare its result to what your code produced. What are the two variables that ``quad`` gives you by default?\n",
    "7. Explain how the trapezoidal rule and Simpson's rule are different. How do the compute the integrals? Explain how the [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature), which is what ``quad`` and its cousins in the ``scipy.integrate`` library do, is different from both of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}