{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Scientific Computing\n", "### Lecture 13: Numerical Integration\n", "#### J.R. Gladden, University of Mississippi, Dept. of Physics" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Integration is the process of computing the \"area bounded by the function and the abscissa (horizontal axis). It can be '+' or '-'.\n", "\n", "Newton-Cotes methods approximate this area by breaking the space up into regular geometric shapes (typically called panels) and computing and adding up the area of all the panels.\n", "\n", "### Trapezoid Method\n", "A panel has the shape of a trapezoid - a rectangle with a triangle on top. Bounding points are connected by a straight line. The error scales like $ \\sim \\mathcal{O}(h^2) $ where $h$ is the width of the rectangle.\n", "#### Area of a trapezoid: $ A = [f(a+h)+f(a)] \\frac{h}{2} $" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import time as t \n", "\n", "def trap(f,xmin,xmax,numPanels):\n", " Isum=0.\n", " h=(xmax-xmin)/float(numPanels)\n", " for i in range(numPanels):\n", " x=xmin+i*h\n", " Isum+=(f(x+h)+f(x))*h/2.\n", " return Isum" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Simpson's 1/3 Method\n", "Here the top of the panel takes the shape of a parabola, which requires an evaluation at 3 points rather than 2. Error scales like $\\sim \\mathcal{O}(h^4)$ but it is not self starting if evaluated on a grid.\n", "### Integral approximation\n", "$$ \\int_a^b f(x) dx \\simeq \\left[ f(a) + 4f(\\frac{a+b}{2})+f(b) \\right] \\frac{h}{3} $$" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def simp13(f,xmin,xmax,numPanels):\n", " Isum=0.\n", " h=(xmax-xmin)/float(numPanels)\n", " for i in range(0,numPanels,2): # NOTE: skipping every other panel\n", " x=xmin+i*h\n", " Isum+=(f(x)+4*f(x+h)+f(x+2*h))*h/3\n", " return Isum" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Usage and comparison" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==================================================\n", "True value = 2.0000\n", "==================================================\n", "Integral by trapeziod with 6 panels = 1.9540972\n", "Error for Trap is 2.2951 percent\n", "==================================================\n", "Integral by Simpson 1/3 with 6 panels = 2.0008632\n", "Error for Simp is 0.0432 percent\n", "==================================================\n" ] } ], "source": [ "%matplotlib wx\n", "\n", "#def f(x): return np.exp(-x) #x**2-0.1*x**3\n", "#Itrue=-(np.exp(-1.0)-np.exp(0.))\n", "\n", "def f(x): return np.sin(x)\n", "Itrue=2.0\n", "\n", "numPanels=6\n", "Itrap=trap(f,0,np.pi,numPanels)\n", "Isimp=simp13(f,0,np.pi,numPanels)\n", "\n", "errortrap=(abs(Itrap-Itrue)/Itrue*100.)\n", "errorsimp=(abs(Isimp-Itrue)/Itrue*100.)\n", "\n", "print('='*50)\n", "print(\"True value = %2.4f\" % Itrue)\n", "print('='*50)\n", "print(\"Integral by trapeziod with %i panels = %3.7f\" % (numPanels,Itrap))\n", "print(\"Error for Trap is %3.4f percent\" % errortrap)\n", "print('='*50)\n", "print(\"Integral by Simpson 1/3 with %i panels = %3.7f\" % (numPanels,Isimp))\n", "print(\"Error for Simp is %3.4f percent\" % errorsimp)\n", "print('='*50)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Visualization\n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "plt.figure()\n", "x=np.linspace(0,1,50)\n", "plt.fill_between(x,np.exp(-x),y2=0,color='yellow')\n", "plt.plot(x,np.exp(-x),'b-',lw=2)\n", "plt.xlim(-0.2,1.2)\n", "plt.xlabel('$x$',fontsize=18)\n", "plt.ylabel('$f(x)$',fontsize=18)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Monte Carlo Methods\n", "Monte Carlo methods work on an entirely different algorithm. Imagine randomly throwing darts at a plot of the function bounded by a rectangular section. Take the ratio of the number of darts falling below the curve to the total number thrown. With enough darts, this will be roughly equal to the ratio of the area under the curve to the toal area of the box.\n", "\n", "This requires a good random number generator and a lot of darts! It is less efficeint that Newton-Cotes for 1D integrals but more efficient for multidimensional integrals." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def MCint(f,xmin,xmax,n,m):\n", " below=0\n", " BoxArea=m*(xmax-xmin)\n", " for i in range(n):\n", " x=np.random.uniform(xmin,xmax) #generates a single random float number\n", " y=np.random.uniform(0,m)\n", " if y < f(x):\n", " below+=1\n", " area = below/float(n)*BoxArea\n", " return area " ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Vectorization: \n", "This pushes the looping down to the underlying C code and is much more efficient." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def MCint_vec(f,xmin,xmax,n,m):\n", " BoxArea=m*(xmax-xmin)\n", " x=np.random.uniform(xmin,xmax,n) #here generates a vector of random floats and length n\n", " y=np.random.uniform(0,m,n)\n", " below=y[y < f(x)].size #picks out indices of y values less than f(x)\n", " area=below/float(n) * BoxArea\n", " return area" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Visualization \n", "Note the use of conditional slicing!" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "def MCint_plot(f,xmin,xmax,n,m):\n", " x=np.random.uniform(xmin,xmax,n)\n", " y=np.random.uniform(0,m,n)\n", " xf=np.linspace(xmin,xmax,100)\n", " #Split random arrays into under and over function - cool use of conditional slicing!!\n", " underY = y[y=f(x)]\n", " overX = x[y>=f(x)]\n", " \n", " plt.plot(overX,overY,'ko',mfc='white',lw=2)\n", " plt.plot(underX,underY,'ko',mfc='red',lw=2)\n", " plt.plot(xf,f(xf),'b-',lw=3)\n", " plt.ylim(0,m)\n", " plt.xlabel('$x$',fontsize=18)\n", " plt.ylabel('$f(x)$',fontsize=18)\n", " return 0" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Usage Examples" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Itrue = 3.1416\n", "==================================================\n", "Monte Carlo Test\n", "Time for 1e+06 points using scalar method: 1.779e+00 seconds\n", "Time for 1e+06 points using vector method: 6.482e-02 seconds\n", "Vector speedup factor: 27.446\n", "==================================================\n", "Error Scaling\n", "\n", "N \t \t Iest \t \t Error(%)\n", "==================================================\n", " 100 \t 3.204425 \t 2.0000\n", " 1000 \t 3.229557 \t 2.8000\n", " 10000 \t 3.131540 \t 0.3200\n", " 100000 \t 3.130660 \t 0.3480\n", " 1000000 \t 3.144841 \t 0.1034\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/joshgladden/Library/Enthought/Canopy/edm/envs/User/lib/python2.7/site-packages/ipykernel/__main__.py:2: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n", " from ipykernel import kernelapp as app\n", "/Users/joshgladden/Library/Enthought/Canopy/edm/envs/User/lib/python2.7/site-packages/ipykernel/__main__.py:3: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n", " app.launch_new_instance()\n" ] } ], "source": [ "def f(x):\n", " return x*np.sin(x)\n", "\n", "plt.figure()\n", "Itrue=np.pi\n", "\n", "print('Itrue = %2.4f'%Itrue)\n", "\n", "\n", "nvals = [100,1000,int(1e4),int(1e5),int(1e6)]\n", "errs=[]\n", "\n", "# Test timing\n", "N=10**6\n", "startTime = t.time()\n", "Iscalar = MCint(f,0,np.pi,N,2)\n", "scalarTime = t.time() - startTime\n", "startTime = t.time()\n", "Ivector = MCint_vec(f,0,np.pi,N,2)\n", "vectorTime = t.time() - startTime\n", "\n", "print(\"=\"*50)\n", "print(\"Monte Carlo Test\")\n", "print(\"Time for %2.0e points using scalar method: %2.3e seconds\" % (N,scalarTime))\n", "print(\"Time for %2.0e points using vector method: %2.3e seconds\" % (N,vectorTime))\n", "print(\"Vector speedup factor: %2.3f\" % (scalarTime/vectorTime) )\n", "print(\"=\"*50)\n", "print(\"Error Scaling\")\n", "print \n", "print('N \\t \\t Iest \\t \\t Error(%)')\n", "print('='*50)\n", "\n", "#MCint_plot(f,0,np.pi,1000,2)\n", "for n in nvals: \n", " #Iest=MCint(f,0,np.pi,n,2)\n", " Iest=MCint_vec(f,0,np.pi,n,2)\n", " errI = abs(Iest-Itrue)/float(Itrue)*100\n", " print('%8i \\t %2.6f \\t %2.4f' % (n,Iest,errI))\n", " errs.append(errI)\n", "\n", "\n", "plt.loglog(nvals,errs,'b-o')\n", "plt.xlabel('Number of points')\n", "plt.ylabel('Error (%)')\n", "\n", "plt.figure()\n", "MCint_plot(f,0,np.pi,1e4,2)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [] } ], "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.13" } }, "nbformat": 4, "nbformat_minor": 0 }