{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import division\n", "from sympy import *\n", "import time\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "def fixfloat(vec):\n", " return list(map(lambda i: '%0.4f' % i, vec))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 1: Analytical Integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the functions\n", "\n", "$f(x) = -4 \\; x \\; (x - 1) \\label{ffunc}$\n", "\n", "and\n", "\n", "$g(x) = \\frac{8}{9} \\cdot (-16\\;x^4 + 40\\;x^3 -35\\;x^2 + 11\\;x)$\n", "\n", "Compute the antiderivatives and evaluate the integrals.\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import division\n", "from sympy import *\n", "x, y = symbols('x y')\n", "l, i = symbols('l, i', integer=True)\n", "g, h = symbols('f g', cls=Function)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-4*x**3/3 + 2*x**2\n", "2/3\n" ] } ], "source": [ "# define symbolic functions : f, intf\n", "# use integrate(func, var)\n", "\n", "f = -4 * x * (x-1)\n", "\n", "intf = integrate(f, x)\n", "\n", "print(intf)\n", "\n", "print(integrate(f,(x,0,1)))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-128*x**5/45 + 80*x**4/9 - 280*x**3/27 + 44*x**2/9\n", "76/135\n" ] } ], "source": [ "# define symbolic functions: g, intg\n", "# use integrate(func, var)\n", "\n", "g = 8 * (-16*x**4 + 40*x**3 - 35*x**2 + 11*x) / 9\n", "\n", "intg = integrate(g,x)\n", "\n", "print(intg)\n", "\n", "print(integrate(g,(x,0,1)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 2: Composite Trapezoidal Rule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Write a function that approximates the integral via the Composite Trapezoidal Rule.\n", "Complete the function template, with $a$, $b$ and $f$ according to the definition above and $n$ being the number of trapezoids used." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def trapezoidalRule(f, a, b, n):\n", " '''\n", " Implements CTR to approximate the integral of a given function. \n", " F ~= T = (b-a) * (f(a)+f(b))/2\n", " '''\n", " integral = 0.0\n", " dh = (b-a)/float(n) # define interval base length\n", " \n", " # for every interval do...\n", " for i in range(n):\n", " l = a + i*dh\n", " r = l + dh\n", " integral += dh * (f(l) + f(r)) * 0.5\n", " \n", " return integral\n", " " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CTR intf = 0.666666000000000\n", "Analytic = 0.6666666666666666\n", "\n", "CTR intg = 0.562961925926402\n", "Analytic = 0.562962962962963\n" ] } ], "source": [ "integ = trapezoidalRule(lambda t: f.subs(x, t), 0.0, 1.0, 1000)\n", "print('CTR intf =',integ)\n", "print('Analytic = ', float(2/3))\n", "\n", "integ = trapezoidalRule(lambda t: g.subs(x, t), 0.0, 1.0, 1000)\n", "print('\\nCTR intg =',integ)\n", "print('Analytic = ', float(76/135))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 3: Composite Simpson Rule" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the same as in Exercise 2 for the Composite Simpson Rule." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def simpsonRule(f, a, b, n):\n", " '''\n", " Implements SR to approximate the integral of given function.\n", " F ~= S = (b-a) * (f(a) + f(b) + 4*f(m)) / 6, where m = (a+b)/2\n", " '''\n", " \n", " integral = 0.0\n", " dh = (b-a)/float(n)\n", " \n", " for i in range(n):\n", " l = a + i*dh\n", " r = l + dh\n", " m = (l + r) * 0.5\n", " \n", " integral += dh * (f(l) + f(r) + 4*f(m)) / 6\n", " \n", " return integral" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SR intf = 0.666666666666667\n", "Analytic = 0.6666666666666666\n", "\n", "SR intg = 0.562962944000001\n", "Analytic = 0.562962962962963\n" ] } ], "source": [ "integ = simpsonRule(lambda t: f.subs(x, t), 0.0, 1.0, 5)\n", "print('SR intf =',integ)\n", "print('Analytic = ', float(2/3))\n", "\n", "integ = simpsonRule(lambda t: g.subs(x, t), 0.0, 1.0, 50)\n", "print('\\nSR intg =',integ)\n", "print('Analytic = ', float(76/135))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 4: Archimedes' Hierarchical Approach" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this exercise we will use Archimedes' approach to approximate the integral.\n", "\n", "Let $\\vec{u} = [u_0,\\dots,u_{n-1}]^T \\in {\\mathbb R}^n, n = 2^l - 1, l \\in {\\mathbb N}$ a vector of function values with\n", "$u_i = f(x_i = \\frac{i + 1}{2^l})$.\n", "\n", "a)\n", "Write a function that transforms a given vector $\\vec{u} \\in {\\mathbb R}^n$ to a similar vector \n", "$\\vec{v} \\in {\\mathbb R}^n$ containing the hierarchical coefficients needed for Archimedes' quadrature approach.\n", "Use the function template below.\n", "\n", "b)\n", "Having computed the vector $\\vec{v}$ with the hierarchical coefficients, implement a function evaluating\n", "the integral.\n", "\n", "c)\n", "Write a function dehierarchize1d similar to hierarchize1d that computes the inverse of the\n", "transformation above.\n", "\n", "NOTE:\n", "We always assume a ZERO boundary case. And the domain is always $[0,1]$." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Task a)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def hierarchize1d(u, maxlv):\n", " ''' \n", " Basic algorithm:\n", " -----------------\n", " for l = maxlevel..1 :\n", " delta = 2**(maxl - l) # detla is indeed the index interval between mid-point and end-point\n", " for j = index of every point of level l:\n", " v[j] = u[j] - 0.5*(u[j-delta] + u[j+delta]) \n", " '''\n", " \n", " v = list(u)\n", " N = len(u)\n", " \n", " for lv in range(maxlv, 1, -1): # no need to include level=1, bc v[j] = u[j] at level 1\n", " delta = 2**(maxlv - lv) # index interval between mid-point and end-point\n", " first = 2**(maxlv - lv) - 1 # first point index of current level\n", " \n", " for j in range(first, N, 2*delta):\n", " #v[j] = u[j] - 0.5 * (u[j-delta] + u[j+delta]) <--- cannot do this directly\n", " # we need to make sure index not out of bound\n", " \n", " v[j] = u[j]\n", " \n", " if (j-delta >= 0):\n", " v[j] -= 0.5 * u[j-delta]\n", " \n", " if (j+delta < N):\n", " v[j] -= 0.5 * u[j+delta]\n", " \n", " return v" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def hierarchize1d_efficient(u, maxl):\n", "\n", " '''\n", " Bitwise operation:\n", " ------------------\n", " x << y : equivalent to (x * 2**y)\n", " x >> y : equivalent to (x / 2**y)\n", " '''\n", "\n", " '''\n", " Improved algorithm:\n", " -------------------\n", " for l = maxlevel-1..2\n", " \n", " We update the grid points of l+1\n", " We get rid of conditionals inside the for loop\n", " \n", " delta_next = 1 << (maxl-l-1)\n", " first_this = (1 << (maxl-l) - 1)\n", " \n", " for j = index of every point of level l:\n", " ....\n", " '''\n", " \n", " v = list(u) # OUTPUT vector\n", " N = len(u)\n", " \n", " for l in range(maxl-1, 0, -1):\n", " \n", " delta_next = 1 << (maxl-l-1) # delta of next level (l+1)\n", " first_this = (1<<(maxl-l)) - 1 # first point index of current level\n", " \n", " for j in range(first_this, N, delta_next<<1):\n", " v[j-delta_next] -= 0.5 * u[j]\n", " v[j+delta_next] -= 0.5 * u[j]\n", " \n", " return v" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = ['0.1250', '0.2500', '0.3750', '0.5000', '0.6250', '0.7500', '0.8750']\n", "u = ['0.4375', '0.7500', '0.9375', '1.0000', '0.9375', '0.7500', '0.4375']\n", "v = ['0.0625', '0.2500', '0.0625', '1.0000', '0.0625', '0.2500', '0.0625']\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEUlJREFUeJzt3X+s3XV9x/Hnm/6Ydxu2i71EoYWyWaodaQIeiInJZNOt\nhYSCP0LahGwszDozdImuGY0ESc0ytYlGMrZRjUFNlFVCmmusu9kYxMzYpbcW2rWkriKuvZ3jwixL\n5lVKfe+Pcy6cXtqec3pOz/d7P30+kqbn+7mfnO+r3889r3zv93vuaWQmkqSyXFR1AEnS4FnuklQg\ny12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpALNr2rHS5YsyeXLl1e1e0mak/bs2fN8Zo52\nmldZuS9fvpyJiYmqdi9Jc1JE/LibeV6WkaQCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQB3L\nPSK+FBHPRcS/n+HrERH3R8ThiNgXEdcOPqYkqRfdnLk/BKw9y9dvBFa0/mwE/q7/WNLw7R57kJ/c\n92Z++YlF/OS+N7N77MGqI0nnrGO5Z+Z3gP85y5RbgK9k0y5gcUS8aVABpWHYPfYgV++5hzcyxUUB\nb2SKq/fcY8FrzhrENffLgCNt20dbY9Kcsez7WxmJl04ZG4mXWPb9rRUlkvoz1BuqEbExIiYiYmJq\namqYu5bO6pI8/ffjJfn8kJNIgzGIcp8ElrVtL22NvUZmbsvMRmY2Rkc7fqiZNDTPxem/H5+LJUNO\nIg3GIMp9DPjD1rtm3g68mJn/NYDnlYbmyLWbmM6Fp4xN50KOXLupokRSfzp+5G9EfB24AVgSEUeB\nTwALADLz74GdwE3AYeBnwB+fr7DS+XLdug+ym+a190vyeZ6LJRx52yauW/fBqqNJ5yQys5IdNxqN\n9PPcJak3EbEnMxud5vkbqpJUIMtdkgpkuUtSgSx3SSqQ5S5JBer4VkipGzv2TrJ1/BDHjk9z6eIR\nNq1Zya3X+CkUVXAtBJa7BmDH3kk2P7qf6RMnAZg8Ps3mR/cDWCpD5lpohpdl1Let44deKZMZ0ydO\nsnX8UEWJLlyuhWZY7urbsePTPY3r/HEtNMNyV98uXTzS07jOH9dCMyx39W3TmpWMLJh3ytjIgnls\nWrOyokQXLtdCM7yhqr7N3KjzHRrVcy00ww8Ok6Q5xA8Ok6QLmOUuSQWy3CWpQJa7JBXIcpekAlnu\nklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtSgSx3SSqQ5S5J\nBbLcJalAlrskFairco+ItRFxKCIOR8Tdp/n65RHxeETsjYh9EXHT4KNKkrrVsdwjYh7wAHAjsArY\nEBGrZk27B9iemdcA64G/HXRQSVL3ujlzvx44nJnPZOZLwMPALbPmJPD61uNFwLHBRZQk9aqbcr8M\nONK2fbQ11u4+4PaIOArsBD58uieKiI0RMRERE1NTU+cQV5LUjUHdUN0APJSZS4GbgK9GxGueOzO3\nZWYjMxujo6MD2rUkabZuyn0SWNa2vbQ11u5OYDtAZn4PeB2wZBABJUm966bcdwMrIuLKiFhI84bp\n2Kw5/wm8CyAi3kqz3L3uIkkV6VjumfkycBcwDjxN810xByJiS0Ssa037GPCBiHgK+DpwR2bm+Qot\nSTq7+d1MysydNG+Uto/d2/b4IPCOwUaTJJ0rf0NVkgpkuUtSgbq6LKPzZ8feSbaOH+LY8WkuXTzC\npjUrufWa2b9GIF14fG30x3Kv0I69k2x+dD/TJ04CMHl8ms2P7gfwm1gXNF8b/fOyTIW2jh965Zt3\nxvSJk2wdP1RRIqkefG30z3Kv0LHj0z2NSxcKXxv9s9wrdOnikZ7GpQuFr43+We4V2rRmJSML5p0y\nNrJgHpvWrKwokVQPvjb65w3VCs3cGPIdAdKpfG30L6r6lIBGo5ETExOV7FuS5qqI2JOZjU7zvCwj\nSQWy3CWpQJa7JBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJU\nIMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIK1FW5R8TaiDgUEYcj4u4z\nzLktIg5GxIGI+NpgY0qSejG/04SImAc8APw+cBTYHRFjmXmwbc4KYDPwjsz8aURccr4CS5I66+bM\n/XrgcGY+k5kvAQ8Dt8ya8wHggcz8KUBmPjfYmJKkXnRT7pcBR9q2j7bG2l0FXBUR342IXRGx9nRP\nFBEbI2IiIiampqbOLbEkqaNB3VCdD6wAbgA2AF+IiMWzJ2XmtsxsZGZjdHR0QLuWJM3WTblPAsva\ntpe2xtodBcYy80Rm/gj4Ac2ylyRVoJty3w2siIgrI2IhsB4YmzVnB82zdiJiCc3LNM8MMKckqQcd\nyz0zXwbuAsaBp4HtmXkgIrZExLrWtHHghYg4CDwObMrMF85XaEnS2UVmVrLjRqORExMTlexbkuaq\niNiTmY1O8/wNVUkqkOUuSQWy3CWpQJa7JBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwl\nqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIK\nZLlLUoEsd0kq0PyqA/Rjx95Jto4f4tjxaS5dPMKmNSu59ZrLqo4lSUC1HTVny33H3kk2P7qf6RMn\nAZg8Ps3mR/cDWPCSKld1R83ZyzJbxw+9ctBmTJ84ydbxQxUlkqRXVd1Rc7bcjx2f7mlckoap6o6a\ns+V+6eKRnsYlaZiq7qg5W+6b1qxkZMG8U8ZGFsxj05qVFSWSpFdV3VFz9obqzA0J3y0jqY6q7qjI\nzM6TItYCnwfmAV/MzE+dYd77gEeA6zJz4mzP2Wg0cmLirFMkSbNExJ7MbHSa1/GyTETMAx4AbgRW\nARsiYtVp5l0M/Dnwb73HlSQNUjfX3K8HDmfmM5n5EvAwcMtp5n0S+DTw8wHmkySdg27K/TLgSNv2\n0dbYKyLiWmBZZn5rgNkkSeeo73fLRMRFwGeBj3Uxd2NETETExNTUVL+7liSdQTflPgksa9te2hqb\ncTFwNfBERDwLvB0Yi4jXXPDPzG2Z2cjMxujo6LmnliSdVTflvhtYERFXRsRCYD0wNvPFzHwxM5dk\n5vLMXA7sAtZ1ereMJOn86VjumfkycBcwDjwNbM/MAxGxJSLWne+AkqTedfVLTJm5E9g5a+zeM8y9\nof9YkqR+zNmPH5AknZnlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtdkgpkuUtSgSx3SSqQ5S5JBbLc\nJalAlrskFchyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBXIcpekAlnuklQgy71q\n+7bD566G+xY3/963vepEUj342ujL/KoDXND2bYdvfgROTDe3XzzS3AZYfVt1uaSq+drom2fuVXps\ny6vfvDNOTDfHpQuZr42+We5VevFob+PShcLXRt8s9yotWtrbuHSh8LXRN8u9Su+6FxaMnDq2YKQ5\nLl3IfG30zXKv0urb4Ob7YdEyIJp/33y/N4wkXxt9i8ysZMeNRiMnJiYq2bckzVURsSczG53meeYu\nSQWy3CWpQF2Ve0SsjYhDEXE4Iu4+zdc/GhEHI2JfRDwWEVcMPqokqVsdyz0i5gEPADcCq4ANEbFq\n1rS9QCMzVwOPAJ8ZdFBJUve6OXO/Hjicmc9k5kvAw8At7RMy8/HM/Flrcxfgm1ElqULdlPtlwJG2\n7aOtsTO5E/h2P6EkSf0Z6AeHRcTtQAN45xm+vhHYCHD55ZcPcteSpDbdnLlPAsvatpe2xk4REe8G\nPg6sy8xfnO6JMnNbZjYyszE6OnoueSVJXeim3HcDKyLiyohYCKwHxtonRMQ1wIM0i/25wceUJPWi\nY7ln5svAXcA48DSwPTMPRMSWiFjXmrYV+HXgGxHxZESMneHpJElD0NU198zcCeycNXZv2+N3DziX\nJKkP/oaqJBXIcpekAlnuklQgy12SCmS5S1KBLHdJKpDlLkkFstwlqUCWuyQVyHKXpAJZ7pJUIMtd\nkgpkuUtSgSx3SSqQ5S5JBbLcJalAlrskFchyl6QCWe6SVCDLXZIKNLfLfd92+NzVcN/i5t/7tled\n6MLlWtSHa1EfFa7F/KHtadD2bYdvfgROTDe3XzzS3AZYfVt1uS5ErkV9uBb1UfFazN0z98e2vHrQ\nZpyYbo5ruFyL+nAt6qPitZi75f7i0d7Gdf64FvXhWtRHxWsxd8t90dLexnX+uBb14VrUR8VrMXfL\n/V33woKRU8cWjDTHNVyuRX24FvVR8VrM3XJffRvcfD8sWgZE8++b7/emURVci/pwLeqj4rWIzBzK\njmZrNBo5MTFRyb4laa6KiD2Z2eg0b+6euUuSzshyl6QCWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWq\n7H3uETEF/LjDtCXA80OIc67qnK/O2aDe+eqcDczXjzpng+7yXZGZo52eqLJy70ZETHTzZv2q1Dlf\nnbNBvfPVORuYrx91zgaDzedlGUkqkOUuSQWqe7lvqzpAB3XOV+dsUO98dc4G5utHnbPBAPPV+pq7\nJOnc1P3MXZJ0DmpR7hGxNiIORcThiLj7NF//nYj4fkS8HBHvr1m2j0bEwYjYFxGPRcQVNcv3pxGx\nPyKejIh/jYhVdcrXNu99EZERMbR3MnRx7O6IiKnWsXsyIv5kWNm6ydeac1vr++9ARHytLtki4nNt\nx+0HEXF8WNm6zHd5RDweEXtbr92bapbvilaf7IuIJyKi9/++KTMr/QPMA34I/CawEHgKWDVrznJg\nNfAV4P01y/a7wK+2Hn8I+Iea5Xt92+N1wD/WKV9r3sXAd4BdQKMu2YA7gL8Z1vE6h3wrgL3Ab7S2\nL6lLtlnzPwx8qWbHbhvwodbjVcCzNcv3DeCPWo9/D/hqr/upw5n79cDhzHwmM18CHgZuaZ+Qmc9m\n5j7glzXM9nhm/qy1uQsY5n9W2U2+/23b/DVgmDdZOuZr+STwaeDnNcxWlW7yfQB4IDN/CpCZz9Uo\nW7sNwNeHkqypm3wJvL71eBFwrGb5VgH/0nr8+Gm+3lEdyv0y4Ejb9tHWWB30mu1O4NvnNdGpusoX\nEX8WET8EPgN8ZEjZoIt8EXEtsCwzvzXEXND92r6v9aPxIxGxbDjRgO7yXQVcFRHfjYhdEbG2RtmA\n5uUF4EpeLaph6CbffcDtEXEU2Enzp4th6SbfU8B7W4/fA1wcEW/oZSd1KPciRMTtQAPYWnWW2TLz\ngcz8LeAvgXuqzjMjIi4CPgt8rOosZ/BNYHlmrgb+CfhyxXlmm0/z0swNNM+OvxARiytN9FrrgUcy\n82TVQWbZADyUmUuBm4Cvtr4f6+IvgHdGxF7gncAk0NMxrMM/ZhJoPyNa2hqrg66yRcS7gY8D6zLz\nF0PKBr0fu4eBW89rolN1yncxcDXwREQ8C7wdGBvSTdWOxy4zX2hbzy8CbxtCrhndrO1RYCwzT2Tm\nj4Af0Cz7OmSbsZ7hXpKB7vLdCWwHyMzvAa+j+bkuw9DN996xzHxvZl5Ds1vIzN5uSg/rJsJZbi7M\nB56h+aPbzM2F3z7D3IcY7g3VjtmAa2jeHFlRx2PXngu4GZioU75Z859geDdUuzl2b2p7/B5gV52O\nHbAW+HLr8RKaP+q/oQ7ZWvPeAjxL6/dpanbsvg3c0Xr8VprX3IeSs8t8S4CLWo//CtjS836GedDP\n8o+9ieZZxw+Bj7fGttA8Ewa4juZZyv8BLwAHapTtn4H/Bp5s/Rmr2bH7PHCgle3xs5VrFflmzR1a\nuXd57P66deyeah27t9Tp2AFB87LWQWA/sL4u2Vrb9wGfGuYx6+HYrQK+21rbJ4E/qFm+9wP/0Zrz\nReBXet2Hv6EqSQWqwzV3SdKAWe6SVCDLXZIKZLlLUoEsd0kqkOUuSQWy3CWpQJa7JBXo/wEvtuop\n4RunzgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "level = 3\n", "xvec = []\n", "for i in range(1, 2**level):\n", " xvec.append(i * 2**(-level))\n", "\n", "uvec = list(map(lambda t: f.subs(x,t), xvec))\n", "\n", "print(\"x =\",fixfloat(xvec))\n", "print(\"u =\",fixfloat(uvec))\n", "\n", "plt.scatter(xvec, uvec)\n", "\n", "vvec = hierarchize1d(uvec, level)\n", "print(\"v =\",fixfloat(vvec))\n", "\n", "plt.scatter(xvec, vvec)\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Basic implementation v = ['0.0625', '0.2500', '0.0625', '1.0000', '0.0625', '0.2500', '0.0625']\n", "Exec time = 0.0005839999999999179\n", "\n", "Efficient implementation v = ['0.0625', '0.2500', '0.0625', '1.0000', '0.0625', '0.2500', '0.0625']\n", "Exec time = 0.0007289999999997576\n" ] } ], "source": [ "tic = time.clock()\n", "v1 = hierarchize1d(uvec, level)\n", "toc1 = time.clock() - tic\n", "\n", "tic = time.clock()\n", "v2 = hierarchize1d_efficient(uvec, level)\n", "toc2 = time.clock() - tic\n", "\n", "print('Basic implementation v =', fixfloat(v1))\n", "print('Exec time = ', toc1)\n", "\n", "print('\\nEfficient implementation v =', fixfloat(v2))\n", "print('Exec time = ', toc2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Task b)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def computeIntegral(v, maxlv): \n", " integ = 0.0\n", " N = len(v)\n", " \n", " for lv in range(maxlv, 0, -1): # for all levels, including level=1\n", " delta = 2**(maxlv - lv) # index interval between mid-point and end-point\n", " first = 2**(maxlv - lv) - 1 # first point index of current level\n", " half_base = 2**(-lv) # length of half base\n", " \n", " for j in range(first, N, 2*delta):\n", " integ += v[j] * half_base\n", " \n", " return integ" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Archimedes intf = 0.656250000000000\n", "Analytic = 0.6666666666666666\n" ] } ], "source": [ "integ = computeIntegral(vvec, level)\n", "print('Archimedes intf =',integ)\n", "print('Analytic = ', float(2/3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Task c)" ] }, { "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.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }