{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework 1\n", "\n", "Note: Since this is the first homework, detailed template (imports, function declarations, comments etc) are provided. Hence, if you are new to python, you are highly encouaged to play with things (for instance inbuilt functions used in the template) beyond what is asked in the homework.
\n", "
\n", "The ellipsis (...) are provided where you are expected to write your solution but feel free to change the template (not over much) in case this style is not to your taste.
\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Link Okpy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from client.api.notebook import Notebook\n", "ok = Notebook('hw1.ok')\n", "_ = ok.auth(inline = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "#Import in-built functions for different integration techniques\n", "#For reference: https://docs.scipy.org/doc/scipy/reference/integrate.html\n", "from scipy.integrate import quad, fixed_quad, romberg, dblquad\n", "#For plotting\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "\n", "def gaussxw(N):\n", " '''Calculate 'N' position and weights for Gaussian quadrature integration\n", " Returns a tuple of 2 arrays, the first array is the position of points and second\n", " array is the corresponding weights.\n", " '''\n", " a = np.linspace(3, 4*N -1, N)/(4*N+2)\n", " x = np.cos(np.pi*a + 1/(8*N*N*np.tan(a)))\n", " eps = 1e-15\n", " delta = 1.\n", " #calc roots\n", " while delta>eps:\n", " p0 = np.ones(N)\n", " p1 = np.copy(x)\n", " for k in range(1, N):\n", " p0, p1 = p1, ((2*k +1)*x*p1 -k*p0)/(k+1)\n", " dp = (N+1)*(p0 -x*p1)/(1-x**2)\n", " dx = p1/dp\n", " x -= dx\n", " delta = max(abs(dx))\n", " \n", " #calc weights\n", " w = 2*(N+1)**2/(N*N*(1 - x**2)*dp**2)\n", " \n", " return x, w\n", "\n", "def gaussxwab(N, a, b):\n", " '''Calculate 'N' position and weights for Gaussian quadrature integration\n", " between 'a' and 'b'\n", " Returns a tuple of 2 arrays, the first array is the position of points and second\n", " array is the corresponding weights.\n", " '''\n", " x, w = gaussxw(N)\n", " return 0.5*(b-a)*x + 0.5*(b+a), 0.5*(b-a)*w\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 1 - Harmonic Oscillator\n", "\n", "The total energy of a harmonic oscillator is given by
\n", "$$ E = \\frac{1}{2}m \\left(\\frac{dx}{dt}\\right)^2 + V(x) $$\n", "Assuming that the potential $V(x)$ is symmetric about $x=0$ and the amplitude of the oscillator is $a$. Then the equation for the time period is given by
\n", "$$ T = \\sqrt{8m} \\int_0^a \\frac{dx}{\\sqrt{V(a) - V(x)}} $$\n", "1. Suppose the potential is $V(x) = x^4$ and mass of the particle $m=1$, write a function that calculates the period for a given amplitude. Use it to calculate time period of 'a = 2'. Use this amplitude for next 2 parts.\n", "2. Use inbuilt fixed_quad function or the manual 'gaussxwab' function above to calculate the time period for different values of 'N' (number of integration points). Calculate the error in the integral by estimating the difference for 'N' & '2N'.\n", " - Approximately, at what 'N' is the absolute error less than $10^{-6}$ for 'a = 2'\n", " - Use inbuilt 'quad' function that returns an error estimate and compare your answer for 'a = 2' (quad uses a more advanced integration technique)\n", "3. Use the inbuilt romberg function to use Romberg integration. \n", " - Note that a simplistic usage with romberg(func, 0, a) will probably give error or 'nan'. Why?\n", " - Assume that we can tolerate the uncertainitiy of $10^{-5}$ in the position.\n", " - Show and output of 'keyword' show = True for 'a = 2'. Use this to estimate error for divmax = 10. Show your calculation and compare it with the python warning.\n", " - Change divmax to change the number of divisions. How does the accuracy change on going from 10 to 15 divisions.\n", "4. Use the function to make a graph of the period for amplitude rangeing from a=0 to a=2. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 1\n", "\n", "def V(x):\n", " 'Potential'\n", " \n", " return ...\n", "\n", "def timep(x, a):\n", " 'Define the function that needs to be integrated (integrand) to calculate time period'\n", " \n", " return ...\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 2\n", "a = 2\n", "\n", "tquad = quad(..., args = (a, ))\n", "print('Inbuilt Gaussian Quadrature gives time period = ', tquad[0], ' with error = ', tquad[1])\n", "\n", "n = \n", "tquadn = fixed_quad(..., args = (a, ), n = ...)[0]\n", "tquadn2 = fixed_quad(..., args = (a, ), n = ...)[0]\n", "print('\\nFor n = %d, the time period is %0.3f, with error = %0.3e'%(n, tquadn, abs(tquadn2 - tquadn)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Part 3 \n", "The romberg gives error on simplistic call of the romberg(func, 0, a) where a is the amplitude because <...>" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 3\n", "#The position tolerance is 10**-5\n", "\n", "tromb = romberg(timep, 0, a-10**-5, args = (a,), divmax=10,show=True)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Calculate Error\n", "..." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#When divmax = 15\n", "...\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#The new error is\n", "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 4\n", "#Calculate points here\n", "#Create 'a' values in the given range\n", "avals = np.linspace(1e-1, 2, 100)\n", "tvals = np.zeros_like(avals)\n", "\n", "#signature for the manual function above\n", "#Feel free to replace it with call to inbuilt function instead\n", "N = #Number of GQ points\n", "for foo in range(avals.size):\n", " a = avals[foo]\n", " x, w = gaussxwab(N, 0, a) \n", " t = timep(x, a)\n", " tvals[foo] = ...\n", " \n", " \n", "\n", "#Draw Figure here \n", "fig, ax = plt.subplots(1, 1)\n", "ax.plot(..., ...)\n", "ax.set_xlabel('Amplitude')\n", "ax.set_ylabel('Time Period')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 2 - Black Body Radiation\n", "\n", "The total rate at which energy is radiated by a black body per unit area over all frequencies is \n", "$$ W = \\frac{2 \\pi k_B^4T^4}{c^2 {h}^3} \\int _0^\\infty \\frac{x^3}{e^x -1} dx $$\n", "1. Write a function to to evaluate the integral in this expression. You will need to change the variables to go from an infinite range to a finite range. What is the change of variable and new functional form?\n", "2. According to Stefan's law, the total energy given off by a black-body per unit area per second is given by \n", "$$ W = \\sigma T^4 $$\n", "Use the integral to calculate the value of Stefan Boltzmann constant $\\sigma$. Use 'fixed_quad' function to do the integral. How accurate you think the answer is?\n", "3. Inbuilt 'quad' function can support an infinite range for integration. Write another function to do the integration from 0 to $\\infty$ and compare your answer. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solution\n", "Part 1\n", "\n", "The variable to go from range 0 to $\\infty$ to a finite range of [...] is - \n", "...\n", "\n", "The change of the function under this change of variables is -\n", "..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 1-2\n", "def blackbody_var(z):\n", " 'Blackbody spectrum after change of variables'\n", " \n", " return ...\n", "\n", "#Constants\n", "k = 1.38064852e-23 \n", "h = 6.626e-34\n", "pi= np.pi\n", "c = 3e8\n", "hb = h /2/pi\n", "prefactor = k**4/c**2/hb**3/4/pi**2 \n", "#True value\n", "stfconst = 5.670367e-8\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Part 2\n", "\n", "#Gaussian Quadrature\n", "N = ...\n", "...fixed_quad(...)\n", "...\n", "...\n", "#Calculate value\n", "stefan = ...\n", "...\n", "error ...\n", "\n", "print(\"Evaluated value = \", stefan, \" with error = \",error)\n", "\n", "print(\"Evaluated value = \", stefan, \"\\nTrue value = \",stfconst, \"\\nResidual from truth = \", (stfconst-stefan)/stfconst)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def blackbody(x):\n", " 'Blackbody Spectrum'\n", " return ...\n", "\n", "stefan2 = quad(..., 0, np.inf)\n", "print(...)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem 3 - Gravitational Pull of Uniform Sheet\n", "\n", "The gravitational force due to a plate felt by a point mass of 1 kg a distnace $z$ from the center of the square in the direction perpendicular to the sheet is given by \n", "$$ F_z = G \\sigma z \\int \\int_{-L/2}^{L/2} \\frac{dx dy}{(x^2 + y^2 + z^2)^{3/2}}$$\n", "where $G = 6.674 \\times 10 ^{-11} m^3 kg^{-1} s^{-2}$ and $\\sigma$ is the mass per unit area.
\n", "\n", "Write a program to calcualte and plot the force as a function of $z$ from $z=0$ to $z=10$ for a sheet of 10 metric tonnes and the sheet is $10\\ m$ on side. Use Gaussian quadrature for the double integral. Though there is a 'dblquad' routine in python, we will make use of the manual functions defined - 'gaussxwab' and 'gaussx' defined above. \n", "- Study how the number of integration points 'N' affects the integral here.\n", "- Based on the above, what do you expect the correct curve to look like. Why do you physically expect that to be the case? Try to make connection with what you might have learnt in electromagnetism. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def force(x, y, z):\n", " 'Write the functional form of force here'\n", " return ...\n", " \n", "#Factors\n", "G = ...\n", "M = 1e4\n", "L = 10.\n", "\n", "#points in z direction\n", "zz = np.logspace(-2, 1, 100)\n", "f, f2 = np.zeros_like(zz), np.zeros_like(zz)\n", "\n", "#Number of points for the integral defining\n", "#points in x,y direction\n", "\n", "N1 = ...\n", "xx1, w1 = gaussxwab(N1, ..., ...)\n", "\n", "for foo in range(zz.size):\n", " z = zz[foo]\n", " ...\n", " f[foo] = ...\n", "\n", "N2 = ...\n", "...\n", "\n", "for foo in range(zz.size):\n", " ...\n", " f2[foo] = \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Hint: \n", "The loop is supposed to calculate the double integral.
\n", "To fill in the for loop, study what is returned by the 'gaussxwab' function and how do you use it to evaluate the integral. The easiest thing that can possibly be done is write a triple for loop.
\n", "However for loops are quite slow in python. You should be able to reduce it to a double loop by using inbuilt functions on numpy array. Take inspiration from here https://docs.scipy.org/doc/numpy/reference/routines.math.html
\n", "The next way to avoid loops in python is use [broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) and/or use of numpy.einsum. Though you are not required to use it here, its nevertheless a handy thing to know about and this problem is one of the simplest cases where you can use it. You will not need to add any new loops inside the one given (though you may need to declare new variables), and it should be faster than both 2 and 3 for loops." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#Make plot\n", "fig, ax = plt.subplots(1, 1)\n", "ax.plot(..., label = ...)\n", "ax.plot(...)\n", "ax.plot(...)\n", "ax.legend(loc = 0)\n", "ax.set_xlabel('z')\n", "ax.set_ylabel('Force')\n", "ax.set_xscale('log')\n", "ax.set_title(\"Force due to a sheet\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Answer:\n", "Write your observation about different 'N' values here. What do you expect the correct curve to look like. Why do you physically expect that to be the case? Try to make connection with what you might have learnt in electromagnetism. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## To Submit\n", "Execute the following cell to submit.\n", "If you make changes, execute the cell again to resubmit the final copy of the notebook, they do not get updated automatically.
\n", "__We recommend that all the above cells should be executed (their output visible) in the notebook at the time of submission.__
\n", "Only the final submission before the deadline will be graded. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "_ = ok.submit()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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": 1 }