PHYS 3310 - Math Methods HW #1 (Fall 2015) completed in Python rather than Mathematica.

phys_3310_hw_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PHYS 3310 - Mathematical Methods in the Physical Sciences\n",
    "\n",
    "## Ch. 1: 1.2.5, 1.2.6, 1.2.10, 1.3.2. 1.4.1, 1.5.1\n",
    "\n",
    "\n",
    "\n",
    "### SymPy\n",
    "\n",
    "SymPy is a powerful, yet completely free, symbolic math library for Python. In conjunction with an IPython notebook it makes an exceptional tool for doing math and physics homework.\n",
    "\n",
    "We start by importing the SymPy module and initializing the LaTeX printer. (Not currently working.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy\n",
    "\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.5 \n",
    "<blockquote>\n",
    "Obtain numerical values, precise to 30 decimal digits, for π^2 and sin 0.1π.\n",
    "<br>\n",
    "<i>Note</i>. Obtain these 30-digit results even if your computer does not normally compute at this precision.\n",
    "</blockquote>\n",
    "\n",
    "SymPy has a variety of built-in functions and constants that we can make ample use of."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi**2 = 9.86960440108935861883449099988\n",
      "sin(0.1*pi) = 0.309016994374947440688093732241\n"
     ]
    }
   ],
   "source": [
    "expressions = [sympy.pi ** 2,\n",
    "               sympy.sin(0.1 * sympy.pi)]\n",
    "\n",
    "for expr in expressions:\n",
    "    print(expr, '=', expr.evalf(30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.6\n",
    "<blockquote>\n",
    "(a) Verify that your computer system knows the special function corresponding\n",
    "to the following integral:\n",
    "<br><br>\n",
    "$$erf(x) = \\frac{2}{\\sqrt{\\pi}} \\int_{0}^{x} e^{-t^{2}}dt$$\n",
    "</blockquote>\n",
    "\n",
    "SymPy allows you to write and manipulate expressions, however you must first define your variables with SymPy's Symbol class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, t = sympy.Symbol('x'), sympy.Symbol('t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With multiple symbols this is done even easier with SymPy's symbols factory function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x, t = sympy.symbols('x t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expressions are fairly straightforward. We'll use the Integral class to symbolically represent our integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2*Integral(exp(-t**2), (t, 0, 18.7206051190171))/sqrt(pi)\n"
     ]
    }
   ],
   "source": [
    "erf = 2 / sympy.sqrt(sympy.pi) * sympy.Integral(sympy.exp(-t**2), (t, 0, x))\n",
    "\n",
    "print(erf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see if SymPy recognizes our expression we'll evaluate our integral with the doit() method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "erf.doit() == sympy.erf(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<blockquote>\n",
    "(b) Write symbolic code to evaluate the integral defining erf(x), and check\n",
    "your code by comparing its output with calls to the function in your\n",
    "symbolic algebra system. Check values for x = 0, x = 1, and x = ∞.\n",
    "</blockquote>\n",
    "\n",
    "For this we'll use the subs() method. subs() takes a pair of arguments in the form subs(symbol, substitute), or an iterable of (old, new) pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "erf(0) = 0\n",
      "erf(1) = 0.842700792949715\n",
      "erf(oo) = 1.00000000000000\n"
     ]
    }
   ],
   "source": [
    "nums = [0, 1, sympy.oo]    # float('inf') will lead to the wrong answer\n",
    "\n",
    "format_str = 'erf({}) = {}'\n",
    "\n",
    "for n in nums:\n",
    "    value = erf.subs(x, n).evalf()\n",
    "    print(format_str.format(n, value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2.10\n",
    "<blockquote>\n",
    "(a) Evaluate the indefinite integral\n",
    "<br><br>\n",
    "$$\\int_{}^{} \\frac{dz}{(z^2 + 4)^{7/2}}$$\n",
    "</blockquote>\n",
    "\n",
    "Sometimes SymPy's algorithm will need a hint to evaluate an integral. Try passing the kwarg manual=True to simulate a more by-hand integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16*4**(-3.5)*z*(z**4 + 10*z**2 + 30)/(15*(z**2 + 4)**(5/2))"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = sympy.Symbol('z')\n",
    "\n",
    "integrand = 1 / (z ** 2 + 4) ** (7/2)\n",
    "\n",
    "sympy.integrate(integrand, z, manual=True).factor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<blockquote>\n",
    "(b) Evaluate the above integral numerically for the range (0, π).\n",
    "</blockquote>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00826264068897049"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.integrate(integrand, (z, 0, sympy.pi), manual=True).evalf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.2\n",
    "<blockquote>\n",
    "Write code that examines the values of two numerical quantities x and y, and\n",
    "produces as output one of the strings Case I or Case II, whichever applies.\n",
    "Case I is when x + y > 10 and 0 < xy < 3; Case II occurs otherwise.\n",
    "</blockquote>\n",
    "\n",
    "This is a simple case that can be done in pure Python using a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compare(x, y):\n",
    "    '''Compares both numeric args and returns the resulting case.'''\n",
    "    \n",
    "    if (x + y > 10) and (0 < x * y < 3):\n",
    "        case = 'Case I'\n",
    "    else:\n",
    "        case = 'Case II'\n",
    "    \n",
    "    return case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test with some random numbers to see if it works like it should."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x=18.1426, y=0.7958: Case II\n",
      "x=14.2853, y=0.5505: Case II\n",
      "x=17.8506, y=0.0532: Case I\n",
      "x=10.1725, y=0.6134: Case II\n",
      "x=16.4024, y=0.1837: Case II\n",
      "x=17.1644, y=0.6150: Case II\n",
      "x=11.9795, y=0.6827: Case II\n",
      "x=13.2275, y=0.3541: Case II\n",
      "x=11.0291, y=0.6018: Case II\n",
      "x=18.7206, y=0.8617: Case II\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "rand_x = (random.uniform(10, 20) for i in range(10))\n",
    "rand_y = (random.uniform(0, 1) for i in range(10))\n",
    "\n",
    "format_str = 'x={:.4f}, y={:.4f}: {}'\n",
    "\n",
    "for x, y in zip(rand_x, rand_y):\n",
    "    print(format_str.format(x, y, compare(x, y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.1\n",
    "<blockquote>\n",
    "Enter the code for the procedure dfac of Example 1.4.3 and verify that it\n",
    "gives correct results for all integer n from −1 to 8.\n",
    "</blockquote>\n",
    "\n",
    "We can use a recursive function, a function that returns itself under some condition, to calculate a double factorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1!! = 1\n",
      "0!! = 1\n",
      "1!! = 1\n",
      "2!! = 2\n",
      "3!! = 3\n",
      "4!! = 8\n",
      "5!! = 15\n",
      "6!! = 48\n",
      "7!! = 105\n",
      "8!! = 384\n"
     ]
    }
   ],
   "source": [
    "def dfac(n):\n",
    "    '''Recursively computes the double factorial of n (i.e., n!!).'''\n",
    "    \n",
    "    if n <= 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return n * dfac(n - 2)\n",
    "\n",
    "format_str = '{}!! = {}'\n",
    "    \n",
    "for i in range(-1, 9):\n",
    "    print(format_str.format(i, dfac(i)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5.1\n",
    "<blockquote>\n",
    "Using the procedure for generating a table of function values (Example 1.5.2), make a table giving values of the function erfc(x) for 0 ≤ x ≤ 10 in unit steps. This function is defined as 1 − erf(x) and is a known function in both maple and mathematica.\n",
    "</blockquote>\n",
    "\n",
    "We'll tackle this problem with Python's math module, part of the standard library. SymPy also has an erfc function built in, but I prefer to use the standard library whenever possible.\n",
    "\n",
    "Using a mixture of range and a list comprehension we can build a table of values with ease."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 1.0), (1, 0.157299207050285), (2, 0.004677734981047268), (3, 2.2090496998585482e-05), (4, 1.541725790028002e-08), (5, 1.5374597944280341e-12), (6, 2.1519736712498925e-17), (7, 4.1838256077794183e-23), (8, 1.1224297172982926e-29), (9, 4.137031746513816e-37), (10, 2.0884875837625433e-45)]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "table = [(x, math.erfc(x)) for x in range(0, 11)]\n",
    "\n",
    "print(table)"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}