{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Integration via von Neumann Rejection

\n", "

\n", " The area under a curve is found via Neumann rejection. The fraction of random points that occur below the integrand is proportional to the value of the integral." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/html": [ "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\"\"\"From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\" \n", "by RH Landau, MJ Paez, and CC Bordeianu (deceased). Copyright R Landau, \n", "Oregon State Unv, MJ Paez, Univ Antioquia, C Bordeianu (deceased), \n", "Univ Bucharest, 2020. Please respect copyright & acknowledge our work.\"\"\" \n", "\n", "# vonNrejectVP: Monte-Carlo integration via vonNeumann rejection (stone throwing)\n", "\n", "import random\n", "from vpython import *\n", "import numpy as np\n", "\n", "N = 100 # points to plot the function\n", "graph = canvas(width=500,height=500,title='vonNeumann Rejection Int')\n", "xsinx = curve(x=list(range(0,N)), color=color.yellow, radius=0.5) \n", "pts = label(pos=vec(-60, -60,0), text='points=', box=0) # Labels\n", "pts2 = label(pos=vec(-30, -60,0), box=0)\n", "inside = label(pos=vec(30,-60,0), text='accepted=', box=0)\n", "inside2 = label(pos=vec(60,-60,0), box=0)\n", "arealbl = label(pos=vec(-65,60,0), text='area=', box=0)\n", "arealbl2 = label(pos=vec(-35,60,0), box=0)\n", "areanal = label(pos=vec(30,60,0), text='analytical=', box=0)\n", "zero = label(pos=vec(-85,-48,0), text='0', box=0)\n", "five = label(pos=vec(-85,50,0), text='5', box=0)\n", "twopi = label(pos=vec(90,-48,0), text='2pi',box=0)\n", "\n", "def fx (x): return x*sin(x)*sin(x) # Integrand\n", "\t\t\n", "def plotfunc(): # Plot function\n", " incr = 2.0*pi/N\n", " for i in range(0,N):\n", " xx = i*incr\n", " xsinx.append(pos=vec(80.0/pi*xx-80, 20*fx(xx)-50,0))\n", " box = curve(pos=[vec(-80,-50,0),vec (-80,50,0),\n", " vec(80,50,0),vec(80,-50,0), vec(-80,-50,0)], color=color.white) \n", " \n", "plotfunc() # Box area = h x w =5*2pi\n", "j = 0\n", "Npts = 3001 # Pts inside box\n", "analyt = (pi)**2 # Analytical integral\n", "areanal.text = 'analytical=%8.5f'%analyt \n", "genpts = points(radius=2)\n", "for i in range(1,Npts): # points inside box\n", " rate(500) # slow process\n", " x = 2.0*np.pi*random() \n", " y = 5*random() \n", " xp = x*80.0/np.pi-80\n", " yp = 20.0*y-50 \n", " pts2.text = '%4s' %i\n", " if y <= fx(x): # Below curve\n", " j += 1 \n", " genpts.append(pos=vec(xp,yp,0), color=color.cyan)\n", " inside2.text='%4s'%j \n", " else: genpts.append(pos=vec(xp,yp,0), color=color.green)\n", " boxarea = 2.0*pi*5.0 \n", " area = boxarea*j/(Npts-1) \n", " arealbl2.text = '%8.5f'%area " ] } ], "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": 2 }