{
"metadata": {
"name": "Labs17_19_sympy"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lab 17\n",
"\n",
"Key by Ryan Brunt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 1"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"\"\"\"\n",
"lab17key.py - written by Ryan Brunt\n",
"\"\"\"\n",
"import numpy as np\n",
"\n",
"\n",
"def deriv(order, accuracy, type, func, vec, h=2e-5):\n",
" \"\"\"\n",
" Return the numerical derivative of a function\n",
"\n",
" Parameters\n",
" ----------\n",
" order : integer\n",
" order of derivative. must be 1 or 2\n",
"\n",
" accuracy : integer\n",
" number of terms in approximations. must be 1,2,3 for forward,\n",
" or 2,4,6 for centered\n",
"\n",
" type : string\n",
" type of numerical derivative to be calculated. must be 'forward'\n",
" or 'centered'\n",
"\n",
" func : function\n",
" function whose derivative you wish to calculate\n",
"\n",
" vec : array_like, dtype=float\n",
" points to evaluate the derivative at\n",
"\n",
" h : float, optional(default=2e-5)\n",
" Step size\n",
"\n",
" Returns\n",
" -------\n",
" a : array_like\n",
" vector of derivative values calculated at the points in `vec`\n",
"\n",
" \"\"\"\n",
" try:\n",
" a = []\n",
" q = vec.shape[0]\n",
" if type == 'forward':\n",
" if order == 1:\n",
" if accuracy == 1:\n",
" for i in range(q):\n",
" p = ((-func(vec[i])+func(vec[i]+h))/h)\n",
" a.append(p)\n",
" elif accuracy == 2:\n",
" for i in range(q):\n",
" p = (((-3./2)*func(vec[i])+2*func(vec[i]+h)-(1./2)*func(vec[i]+2*h))/(h))\n",
" a.append(p)\n",
" elif accuracy ==3:\n",
" for i in range(q):\n",
" p = (((-11./6)*func(vec[i])+3*func(vec[i]+h)-(3./2)*func(vec[i]+2*h))+(1./3)*func(vec[i]+3*h))/h\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy: must be 1, 2, or 3'\n",
" elif order == 2:\n",
" if accuracy == 1:\n",
" for i in range(q):\n",
" p = (func(vec[i])-2*func(vec[i]+h)+func(vec[i]+2*h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==2:\n",
" for i in range(q):\n",
" p = (2*func(vec[i])-5*func(vec[i]+h)+4*func(vec[i]+2*h)-func(vec[i]+3*h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==3:\n",
" for i in range(q):\n",
" p = ((35./12)*func(vec[i])-(26./3)*func(vec[i]+h)+(19./2)*func(vec[i]+2*h)-(14./3)*func(vec[i]+3*h)+(11./12)*func(vec[i]+4*h))/(h**2)\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy: must be 1 2 or 3'\n",
" else:\n",
" print 'invalid order'\n",
" elif type == 'centered':\n",
" if order == 1:\n",
" if accuracy == 2:\n",
" for i in range(q):\n",
" p = ((-1./2)*func(vec[i]-h)+(1./2)*func(vec[i]+h))/h\n",
" a.append(p)\n",
" elif accuracy ==4:\n",
" for i in range(q):\n",
" p = ((1./12)*func(vec[i]-2*h)-(2./3)*func(vec[i]-h)+(2./3)*func(vec[i]+h)-(1./12)*func(vec[i]+2*h))/h\n",
" a.append(p)\n",
" elif accuracy ==6:\n",
" for i in range(q):\n",
" p = ((-1./60)*func(vec[i]-3*h)+(3./20)*func(vec[i]-2*h)-(3./4)*func(vec[i]-h)+(3./4)*func(vec[i]+h)-(3./20)*func(vec[i]+2*h)+(1./60)*func(vec[i]+3*h))/h\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy:must be 2,4, or 6'\n",
" elif order == 2:\n",
" if accuracy == 2:\n",
" for i in range(q):\n",
" p = (func(vec[i]-h)-2*func(vec[i])+func(vec[i]+h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==4:\n",
" for i in range(q):\n",
" p = ((-1./12)*func(vec[i]-2*h)+(4./3)*func(vec[i]-h)-(5./2)*func(vec[i])+(4./3)*func(vec[i]+h)-(1./12)*func(vec[i]+2*h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==6:\n",
" for i in range(q):\n",
" p = ((1./90)*func(vec[i]-3*h)-(3./20)*func(vec[i]-2*h)+(3./2)*func(vec[i]-h)-(49./18)*func(vec[i])+(3./2)*func(vec[i]+h)-(3./20)*func(vec[i]+2*h)+(1./90)*func(vec[i]+3*h))/(h**2)\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy: must be 2, 4, or 6'\n",
" else:\n",
" print 'invalid order'\n",
" elif type == 'backward':\n",
" if order == 1:\n",
" if accuracy == 1:\n",
" for i in range(q):\n",
" p = (func(vec[i])-func(vec[i]-h))/h\n",
" a.append(p)\n",
" elif accuracy ==2:\n",
" for i in range(q):\n",
" p = ((3./2)*func(vec[i])-2*func(vec[i]-h)+(1./2)*func(vec[i]-2*h))/(h)\n",
" a.append(p)\n",
" elif accuracy ==3:\n",
" for i in range(q):\n",
" p = (((11./6)*func(vec[i])-3*func(vec[i]-h)+(3./2)*func(vec[i]-2*h))-(1./3)*func(vec[i]-3*h))/h\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy: must be 1, 2, or 3'\n",
" elif order == 2:\n",
" if accuracy == 1:\n",
" for i in range(q):\n",
" p = (func(vec[i])-2*func(vec[i]-h)+func(vec[i]-2*h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==2:\n",
" for i in range(q):\n",
" p = (2*func(vec[i])-5*func(vec[i]-h)+4*func(vec[i]-2*h)-func(vec[i]-3*h))/(h**2)\n",
" a.append(p)\n",
" elif accuracy ==3:\n",
" for i in range(q):\n",
" p = ((35./12)*func(vec[i])-(26./3)*func(vec[i]-h)+(19./2)*func(vec[i]-2*h)-(14./3)*func(vec[i]-3*h)+(11./12)*func(vec[i]-4*h))/(h**2)\n",
" a.append(p)\n",
" else:\n",
" print 'invalid accuracy: must be 1, 2, or 3'\n",
" else:\n",
" print 'invalid order'\n",
" else:\n",
" print 'invalid type'\n",
" return np.asarray(a)\n",
" except AttributeError:\n",
" print 'you n00b. you didnt enter all the parameters'\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 2"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import pandas as pd\n",
"from itertools import product\n",
"\n",
"\n",
"def fun17_2(x):\n",
" \"\"\"\n",
" Analytical derivative is -3 * sin(x)\n",
" \"\"\"\n",
" return 3 * np.cos(x)\n",
"\n",
"\n",
"def dfunc17_2(x, order=1):\n",
" \"\"\"derivative of fun17_2\"\"\"\n",
" if order == 1:\n",
" return -3 * np.sin(x)\n",
" elif order == 2:\n",
" return -3 * np.cos(x)\n",
"\n",
"test_x = np.arange(-2, 2, .025)\n",
"\n",
"orders = [1, 2]\n",
"acc = [1, 2, 3]\n",
"h_vals = [1, 1e-3, 1e-5, 1e-7, 1e-9]\n",
"der_types = pd.MultiIndex.from_tuples(list(product(orders, acc, h_vals)))\n",
"results = pd.DataFrame(index=der_types, columns=['centered', 'forward'])\n",
"results.columns.names = ['errors_by_type']\n",
"results.index.names = ['Order', 'Accuracy', 'Step (h)']\n",
"\n",
"for ord_i in orders:\n",
" actual = dfunc17_2(test_x, ord_i)\n",
" for acc_i in acc:\n",
" for h_i in h_vals:\n",
" cent = deriv(ord_i, 2 * acc_i, 'centered', fun17_2, test_x, h=h_i)\n",
" forward = deriv(ord_i, acc_i, 'forward', fun17_2, test_x, h=h_i)\n",
" \n",
" results.ix[ord_i, acc_i, h_i]['centered'] = np.abs(actual - cent).max()\n",
" results.ix[ord_i, acc_i, h_i]['forward'] = np.abs(actual - forward).max()\n",
"\n",
"results"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"
\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" errors_by_type | \n",
" centered | \n",
" forward | \n",
"
\n",
" \n",
" Order | \n",
" Accuracy | \n",
" Step (h) | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" 1 | \n",
" 1 | \n",
" 1 | \n",
" 0.4755828 | \n",
" 1.458758 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 4.999959e-07 | \n",
" 0.0015 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 7.426371e-11 | \n",
" 1.5e-05 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 3.758555e-09 | \n",
" 1.525071e-07 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 3.465643e-07 | \n",
" 6.126496e-07 | \n",
"
\n",
" \n",
" 2 | \n",
" 1 | \n",
" 0.08876399 | \n",
" 0.9332589 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 6.346035e-13 | \n",
" 9.999933e-07 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 5.761036e-11 | \n",
" 1.885256e-10 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 4.868778e-09 | \n",
" 1.073327e-08 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 5.516318e-07 | \n",
" 1.320662e-06 | \n",
"
\n",
" \n",
" 3 | \n",
" 1 | \n",
" 0.0176361 | \n",
" 0.6691347 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 7.482903e-13 | \n",
" 7.512357e-10 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 7.683065e-11 | \n",
" 2.793474e-10 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 7.228002e-09 | \n",
" 1.948056e-08 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 6.662133e-07 | \n",
" 2.166962e-06 | \n",
"
\n",
" \n",
" 2 | \n",
" 1 | \n",
" 1 | \n",
" 0.2418138 | \n",
" 2.768717 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 2.503239e-07 | \n",
" 0.002999979 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 9.967765e-06 | \n",
" 3.802239e-05 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 0.1004527 | \n",
" 0.106455 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 890.8112 | \n",
" 1333.889 | \n",
"
\n",
" \n",
" 2 | \n",
" 1 | \n",
" 0.03049187 | \n",
" 2.413696 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 1.609021e-09 | \n",
" 2.750664e-06 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 1.758813e-05 | \n",
" 3.128172e-05 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 0.1721033 | \n",
" 0.3668233 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 2107.803 | \n",
" 3776.379 | \n",
"
\n",
" \n",
" 3 | \n",
" 1 | \n",
" 0.004586739 | \n",
" 2.09567 | \n",
"
\n",
" \n",
" 0.001 | \n",
" 2.134336e-09 | \n",
" 9.640946e-09 | \n",
"
\n",
" \n",
" 1e-05 | \n",
" 2.682928e-05 | \n",
" 9.166569e-05 | \n",
"
\n",
" \n",
" 1e-07 | \n",
" 0.2289471 | \n",
" 0.7884678 | \n",
"
\n",
" \n",
" 1e-09 | \n",
" 2158.798 | \n",
" 8217.208 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"output_type": "pyout",
"prompt_number": 3,
"text": [
"errors_by_type centered forward\n",
"Order Accuracy Step (h) \n",
"1 1 1 0.4755828 1.458758\n",
" 0.001 4.999959e-07 0.0015\n",
" 1e-05 7.426371e-11 1.5e-05\n",
" 1e-07 3.758555e-09 1.525071e-07\n",
" 1e-09 3.465643e-07 6.126496e-07\n",
" 2 1 0.08876399 0.9332589\n",
" 0.001 6.346035e-13 9.999933e-07\n",
" 1e-05 5.761036e-11 1.885256e-10\n",
" 1e-07 4.868778e-09 1.073327e-08\n",
" 1e-09 5.516318e-07 1.320662e-06\n",
" 3 1 0.0176361 0.6691347\n",
" 0.001 7.482903e-13 7.512357e-10\n",
" 1e-05 7.683065e-11 2.793474e-10\n",
" 1e-07 7.228002e-09 1.948056e-08\n",
" 1e-09 6.662133e-07 2.166962e-06\n",
"2 1 1 0.2418138 2.768717\n",
" 0.001 2.503239e-07 0.002999979\n",
" 1e-05 9.967765e-06 3.802239e-05\n",
" 1e-07 0.1004527 0.106455\n",
" 1e-09 890.8112 1333.889\n",
" 2 1 0.03049187 2.413696\n",
" 0.001 1.609021e-09 2.750664e-06\n",
" 1e-05 1.758813e-05 3.128172e-05\n",
" 1e-07 0.1721033 0.3668233\n",
" 1e-09 2107.803 3776.379\n",
" 3 1 0.004586739 2.09567\n",
" 0.001 2.134336e-09 9.640946e-09\n",
" 1e-05 2.682928e-05 9.166569e-05\n",
" 1e-07 0.2289471 0.7884678\n",
" 1e-09 2158.798 8217.208"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lab 19\n",
"\n",
"Key by Spencer Lyon"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from itertools import product\n",
"from time import time\n",
"import numpy as np\n",
"import pandas as pd"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 1"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def jacobian(f, x0, h=1e-8, how='centered', *args, **kwargs):\n",
" \"\"\"\n",
" Numerically estimate the Jacobian matrix for the function f.\n",
" If f operates :math:`\\\\mathbb{R}^n \\\\rightarrow \\\\mathbb{R}^m`, the\n",
" return value of this function will be an (n, m) matrix of\n",
" partial derivatives.\n",
"\n",
" Parameters\n",
" ----------\n",
" f : function\n",
" The function for which you would like to calculate the jacobian.\n",
" Must accept as an input a numpy array of the same size as\n",
" x0 (also can accept other args passed as *args or **kwargs) and\n",
" return a numpy array of values.\n",
"\n",
" x0 : array_like, dtype=float, shape=(n,)\n",
" The array of values representing the point where the jacobian\n",
" is to be evaluated.\n",
"\n",
" h : float, optional(default=1e-8)\n",
" The step size to be used in calculating the derivative.\n",
"\n",
" how : str\n",
" A string specifying the method of differentiation. Acceptable\n",
" values are centered (c), forwards (f), backward (b).\n",
"\n",
" *args, **kwargs :\n",
" Other arguments that should be passed to f for evaluation.\n",
"\n",
" Returns\n",
" -------\n",
" jac : array_like, dtype=float, shape=(n, m)\n",
" The numpy array of partial derivatives representing the jacobian\n",
"\n",
" \"\"\"\n",
" # Make sure we have an array\n",
" x0 = np.asarray(x0)\n",
"\n",
" # Get n and m\n",
" f0 = f(x0, *args, **kwargs)\n",
" nx = x0.size\n",
" mx = f0.size\n",
"\n",
" jac = np.zeros((nx, mx))\n",
"\n",
" ident = np.eye(nx)\n",
"\n",
" if how.startswith('f'):\n",
" for xi in xrange(nx):\n",
" ei = ident[:, xi]\n",
" jac[xi, :] = (f(x0 + h*ei, *args, **kwargs) - f0) / h\n",
"\n",
" elif how.startswith('c'):\n",
" for xi in xrange(nx):\n",
" ei = ident[:, xi]\n",
" jac[xi, :] = (f(x0 + h*ei, *args, **kwargs) -\n",
" f(x0 - h*ei, *args, **kwargs)) / (2 * h)\n",
"\n",
" elif how.startswith('b'):\n",
" for xi in xrange(nx):\n",
" ei = ident[:, xi]\n",
" jac[xi, :] = (f0 - f(x0 - h*ei, *args, **kwargs)) / h\n",
"\n",
" return jac.T"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def test_fun(x):\n",
" x1, x2 = x\n",
" return np.array([np.exp(x1)*np.sin(x2) + x2 ** 3,\n",
" 3*x2 - np.cos(x1)])\n",
"\n",
"def real_jac(x):\n",
" x1, x2 = x\n",
" return np.array([[np.exp(x1) * np.sin(x2), np.exp(x1)*np.cos(x2) + 3*x2**2],\n",
" [np.sin(x1), 3]])\n",
"\n",
"x_grid = np.linspace(-1, 1, 100)\n",
"y_grid = np.linspace(-1, 1, 100)\n",
"\n",
"max_err = 0.\n",
"\n",
"methods = 'cbf'\n",
"h_vals = [1., 1e-4, 1e-5, 1e-6, 1e-7, 1e-10]\n",
"\n",
"ind = pd.MultiIndex.from_tuples(list(product(list(methods), h_vals)),\n",
" names=['method', 'h'])\n",
"jac_results = pd.DataFrame(index=ind, columns=['max_error', 'time_delta'])\n",
"\n",
"print 'Computing max error in my Jacobian for various how and h values'\n",
"for der_type in product(methods, h_vals):\n",
" how, h = der_type\n",
" max_err = 0.\n",
" for point in product(x_grid, y_grid):\n",
"\n",
" numerical = jacobian(test_fun, point, how=how, h=h)\n",
" analytical = real_jac(point)\n",
" this_err = np.abs(numerical - analytical).max()\n",
" max_err = this_err if this_err > max_err else max_err\n",
" \n",
" # Time methods. Should be same for all points so just use 1\n",
" t1 = time()\n",
" _ = jacobian(test_fun, point, how=how, h=h)\n",
" num_time = time() - t1\n",
" \n",
" t1 = time()\n",
" _ = real_jac(point)\n",
" ana_time = time() - t1\n",
" \n",
" jac_results.ix[how, h]['max_error'] = max_err\n",
" jac_results.ix[how, h]['time_delta'] = num_time - ana_time\n",
"\n",
"jac_results"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Computing max error in my Jacobian for various how and h values\n"
]
},
{
"html": [
"\n",
"
\n",
" \n",
" \n",
" | \n",
" | \n",
" max_error | \n",
" time_delta | \n",
"
\n",
" \n",
" method | \n",
" h | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" c | \n",
" 1.000000e+00 | \n",
" 0.9684898 | \n",
" 0.000177145 | \n",
"
\n",
" \n",
" 1.000000e-04 | \n",
" 9.667494e-09 | \n",
" 0.0001790524 | \n",
"
\n",
" \n",
" 1.000000e-05 | \n",
" 1.175886e-10 | \n",
" 0.0001740456 | \n",
"
\n",
" \n",
" 1.000000e-06 | \n",
" 5.293845e-10 | \n",
" 0.0001809597 | \n",
"
\n",
" \n",
" 1.000000e-07 | \n",
" 6.391792e-09 | \n",
" 0.0001749992 | \n",
"
\n",
" \n",
" 1.000000e-10 | \n",
" 4.380994e-06 | \n",
" 0.0001809597 | \n",
"
\n",
" \n",
" b | \n",
" 1.000000e+00 | \n",
" 3.826186 | \n",
" 0.0001118183 | \n",
"
\n",
" \n",
" 1.000000e-04 | \n",
" 0.0002845317 | \n",
" 0.0001239777 | \n",
"
\n",
" \n",
" 1.000000e-05 | \n",
" 2.845231e-05 | \n",
" 0.0001139641 | \n",
"
\n",
" \n",
" 1.000000e-06 | \n",
" 2.845191e-06 | \n",
" 0.0001130104 | \n",
"
\n",
" \n",
" 1.000000e-07 | \n",
" 2.870152e-07 | \n",
" 0.0001199245 | \n",
"
\n",
" \n",
" 1.000000e-10 | \n",
" 9.935029e-06 | \n",
" 0.0001139641 | \n",
"
\n",
" \n",
" f | \n",
" 1.000000e+00 | \n",
" 3.826186 | \n",
" 0.000109911 | \n",
"
\n",
" \n",
" 1.000000e-04 | \n",
" 0.0002845317 | \n",
" 0.0001111031 | \n",
"
\n",
" \n",
" 1.000000e-05 | \n",
" 2.845231e-05 | \n",
" 0.0001130104 | \n",
"
\n",
" \n",
" 1.000000e-06 | \n",
" 2.845191e-06 | \n",
" 0.0001108646 | \n",
"
\n",
" \n",
" 1.000000e-07 | \n",
" 2.870152e-07 | \n",
" 0.0001130104 | \n",
"
\n",
" \n",
" 1.000000e-10 | \n",
" 8.019226e-06 | \n",
" 0.0001130104 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"output_type": "pyout",
"prompt_number": 3,
"text": [
" max_error time_delta\n",
"method h \n",
"c 1.000000e+00 0.9684898 0.000177145\n",
" 1.000000e-04 9.667494e-09 0.0001790524\n",
" 1.000000e-05 1.175886e-10 0.0001740456\n",
" 1.000000e-06 5.293845e-10 0.0001809597\n",
" 1.000000e-07 6.391792e-09 0.0001749992\n",
" 1.000000e-10 4.380994e-06 0.0001809597\n",
"b 1.000000e+00 3.826186 0.0001118183\n",
" 1.000000e-04 0.0002845317 0.0001239777\n",
" 1.000000e-05 2.845231e-05 0.0001139641\n",
" 1.000000e-06 2.845191e-06 0.0001130104\n",
" 1.000000e-07 2.870152e-07 0.0001199245\n",
" 1.000000e-10 9.935029e-06 0.0001139641\n",
"f 1.000000e+00 3.826186 0.000109911\n",
" 1.000000e-04 0.0002845317 0.0001111031\n",
" 1.000000e-05 2.845231e-05 0.0001130104\n",
" 1.000000e-06 2.845191e-06 0.0001108646\n",
" 1.000000e-07 2.870152e-07 0.0001130104\n",
" 1.000000e-10 8.019226e-06 0.0001130104"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 2"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def hessian(f, x0, h=1e-4, *args, **kwargs):\n",
" \"\"\"\n",
" Numerically estimate the Hessian matrix for the function f.\n",
" It is assumed that f goes from :math:`\\mathbb{R}^n \n",
" \\rightarrow \\mathbb{R}`.\n",
"\n",
" Parameters\n",
" ----------\n",
" f : function\n",
" The function for which you would like to calculate the jacobian.\n",
" Must accept as an input a numpy array of the same size as\n",
" x0 (also can accept other args passed as *args or **kwargs) and\n",
" return a numpy array of values.\n",
"\n",
" x0 : array_like, dtype=float, shape=(n,)\n",
" The array of values representing the point where the jacobian\n",
" is to be evaluated.\n",
"\n",
" h : float, optional(default=1e-4)\n",
" The step size to be used in calculating the derivative.\n",
"\n",
" *args, **kwargs :\n",
" Other arguments that should be passed to f for evaluation.\n",
"\n",
" Returns\n",
" -------\n",
" hes : array_like, dtype=float, shape=(n, n)\n",
" The numpy array of mixed second ordre partial derivatives\n",
" at the point x0\n",
"\n",
" \"\"\"\n",
" # Make sure we have an array\n",
" x0 = np.asarray(x0)\n",
"\n",
" # Get n\n",
" nx = x0.size\n",
"\n",
" hes = np.zeros((nx, nx))\n",
"\n",
" ident = np.eye(nx)\n",
"\n",
" for xi in xrange(nx):\n",
" ei = ident[:, xi]\n",
" for xj in xrange(nx):\n",
" ej = ident[:, xj]\n",
" hes[xi, xj] = (f(x0 + h*(ei + ej), *args, **kwargs) -\n",
" f(x0 + h*(ei - ej), *args, **kwargs) -\n",
" f(x0 + h*(ej - ei), *args, **kwargs) +\n",
" f(x0 - h*(ei + ej), *args, **kwargs)) / (4*h**2)\n",
"\n",
" return hes.T"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def test_hes(vec):\n",
" \"\"\"Rosenbrock banana function\"\"\"\n",
" x, y = vec\n",
" return (1 - x)**2 + 100*(y - x**2)**2\n",
"\n",
"\n",
"def ana_ban(vec):\n",
" x, y = vec\n",
" return np.array([[1200*x**2 - 400*y + 2, -400*x],\n",
" [-400*x, 200]])\n",
"\n",
"h_vals = [1., 1e-4, 1e-5, 1e-6, 1e-7, 1e-10]\n",
"hes_results = pd.DataFrame(index=h_vals, columns=['max_error', 'time_delta'])\n",
"\n",
"hes_x = np.linspace(-2, 2, 100)\n",
"hes_y = np.linspace(0, 2, 100)\n",
"\n",
"print 'Computing max error in my Hessian for various how and h values'\n",
"for h in h_vals:\n",
" max_err = 0.\n",
" for point in product(hes_x, hes_y):\n",
"\n",
" numerical = hessian(test_hes, point, h=h)\n",
" analytical = ana_ban(point)\n",
" this_err = np.abs(numerical - analytical).max()\n",
" max_err = this_err if this_err > max_err else max_err\n",
" \n",
" # Time methods. Should be same for all points so just use 1\n",
" t1 = time()\n",
" _ = hessian(test_hes, point, h=h)\n",
" num_time = time() - t1\n",
" \n",
" t1 = time()\n",
" _ = ana_ban(point)\n",
" ana_time = time() - t1\n",
" \n",
" hes_results.ix[h]['max_error'] = max_err\n",
" hes_results.ix[h]['time_delta'] = num_time - ana_time\n",
"\n",
"hes_results"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Computing max error in my Hessian for various how and h values\n"
]
},
{
"html": [
"\n",
"
\n",
" \n",
" \n",
" | \n",
" max_error | \n",
" time_delta | \n",
"
\n",
" \n",
" \n",
" \n",
" 1.000000e+00 | \n",
" 800 | \n",
" 0.0003399849 | \n",
"
\n",
" \n",
" 1.000000e-04 | \n",
" 3.043648e-05 | \n",
" 0.0003130436 | \n",
"
\n",
" \n",
" 1.000000e-05 | \n",
" 0.0029908 | \n",
" 0.0003321171 | \n",
"
\n",
" \n",
" 1.000000e-06 | \n",
" 0.2581446 | \n",
" 0.0003361702 | \n",
"
\n",
" \n",
" 1.000000e-07 | \n",
" 27.96591 | \n",
" 0.0003869534 | \n",
"
\n",
" \n",
" 1.000000e-10 | \n",
" 2.274182e+07 | \n",
" 0.0003211498 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"output_type": "pyout",
"prompt_number": 5,
"text": [
" max_error time_delta\n",
"1.000000e+00 800 0.0003399849\n",
"1.000000e-04 3.043648e-05 0.0003130436\n",
"1.000000e-05 0.0029908 0.0003321171\n",
"1.000000e-06 0.2581446 0.0003361702\n",
"1.000000e-07 27.96591 0.0003869534\n",
"1.000000e-10 2.274182e+07 0.0003211498"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# sympy lab\n",
"\n",
"Key by Spencer Lyon"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import sympy as sym\n",
"import matplotlib.pyplot as plt\n",
"%load_ext sympyprinting"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 1"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print('7 appears %i times in the first 10000 digits of pi' %\n",
" str(sym.pi.n(10000)).count('7'))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"7 appears 970 times in the first 10000 digits of pi\n"
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 2"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"a, b, c, x, x1, x2, x3, h, y1, y2, y3 = sym.symbols('a, b, c, x, x1, x2, x3, \\\n",
" h, y1, y2, y3')\n",
"x2 = x1 + h\n",
"x3 = x1 + 2 * h\n",
"eqns = (sym.Eq(y1, a + b * x1 + c * x1 ** 2),\n",
" sym.Eq(y2, a + b * x2 + c * x2 ** 2),\n",
" sym.Eq(y3, a + b * x3 + c * x3 ** 2))\n",
"soln = sym.solve(eqns, a, b, c)\n",
"print('(a, b, c) = ')\n",
"soln[a], soln[b], soln[c]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(a, b, c) = \n"
]
},
{
"latex": [
"$$\\begin{pmatrix}\\frac{2 h^{2} y_{1} + 3 h x_{1} y_{1} - 4 h x_{1} y_{2} + h x_{1} y_{3} + x_{1}^{2} y_{1} - 2 x_{1}^{2} y_{2} + x_{1}^{2} y_{3}}{2 h^{2}}, & \\frac{- 3 h y_{1} + 4 h y_{2} - h y_{3} - 2 x_{1} y_{1} + 4 x_{1} y_{2} - 2 x_{1} y_{3}}{2 h^{2}}, & \\frac{y_{1} - 2 y_{2} + y_{3}}{2 h^{2}}\\end{pmatrix}$$"
],
"output_type": "pyout",
"prompt_number": 8,
"text": [
"(\n",
" 2 2 2 2 \n",
"2\u22c5h \u22c5y\u2081 + 3\u22c5h\u22c5x\u2081\u22c5y\u2081 - 4\u22c5h\u22c5x\u2081\u22c5y\u2082 + h\u22c5x\u2081\u22c5y\u2083 + x\u2081 \u22c5y\u2081 - 2\u22c5x\u2081 \u22c5y\u2082 + x\u2081 \u22c5y\u2083\n",
"\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
" 2 \n",
" 2\u22c5h ,\n",
" \n",
"-3\u22c5h\u22c5y\u2081 + 4\u22c5h\u22c5y\u2082 - h\u22c5y\u2083 - 2\u22c5x\u2081\u22c5y\u2081 + 4\u22c5x\u2081\u22c5y\u2082 - 2\u22c5x\u2081\u22c5y\u2083\n",
"\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
" 2 \n",
" 2\u22c5h ,\n",
" \n",
"y\u2081 - 2\u22c5y\u2082 + y\u2083\n",
"\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
" 2 \n",
" 2\u22c5h )"
]
}
],
"prompt_number": 8
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 3"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"y = a + b * x + c * x ** 2\n",
"\n",
"soln[x] = (x1 + x2) / 2\n",
"y_est1 = sym.simplify(y.subs(soln))\n",
"print 'y(x1 + x2) / 2 ->\\n'\n",
"y_est1"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"y(x1 + x2) / 2 ->\n",
"\n"
]
},
{
"latex": [
"$$\\frac{3}{8} y_{1} + \\frac{3}{4} y_{2} - \\frac{1}{8} y_{3}$$"
],
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAHUAAAAeCAYAAAASN1ElAAAABHNCSVQICAgIfAhkiAAAA+5JREFU\naIHt2kmIHFUcx/FP4hJNMgpelMRlMl4UNUJwRU0GNzwEIWEMaFzQiASUJK4H9RAPwhgVRUXBeGgX\nYhBXMBqIhxE8xIMYd4gg8SIGPahxwajEw78aa9rqqeru6upu7S80Xa/eq/r/+q3/9/7NkP8cB6Wu\nlyafs3AbfsJXvRDVIoOkewRbsBM/VmHwe1yXXF+BXxMR/c6g6L4RG3EAo1UZPRXzkusJ/KE/K6eR\nQdNdaaOm2YJ7emG4QwZBd9cb9eCG9BJciF/wSDcNl8yg6q6Um/AB5vfI/mL/7nBF6LXuIlQ2/Z6D\nvViUpE9KjE9UYTyDmmI/vN90F6HrjTo7+f4Tn+GbJD0mHI5d3TReAoOquzKuxu24Fa/g0owyo3gc\n23BVQ97NeKckLTXFe3MR3VSnvRmr8ZQYqVtxS5NyoyrW+SQOxXp83JD3vhBbBjXlT1FVae+UjnXO\nziuQ4ny8h/24DLtTefOEB/puC++rkkHRXrnOYzAHC/EXVqTyLhHTyikl2aopd6RWqb0TStHZyrbh\n2+R7FX7GW6m8C8Rx3ectvA+exekZ948XZ7n7M/LWiG1LK5Sl/TQ8jVkF7e7C2oJlKbmOD+R80mzH\n6w33pvBqw70RvCwaqFVqio3UPN3taD9bOF0bsUMECzqlbJ1LcY04T34BF3cq8As8kErPwW/YkLrX\n6eF1rc3n8sjTPh+TqfxVIkCwsAtaZiJP54xBjFYcpTpf46hUehKHmb6APyMatd/I0z6Gu3Bikt6O\nw3FeVQIT8nSOi1mQGDiHpB9Or6mLsBy/42i8iQ8zDG4QjfaE6CFniNjgR+3/htJ4GK8JDzKLPO2f\n4Fz/xGOPS76/LEFb0fotovPTVNmVYgDty3rRg6nrWSLikcdssbg/1yS/yul3mTgyHC9YPk87PC86\nShm0U78017kEd2Az5jZ7eLeITdbZnFHmRdM3xCuEh3pyk3e226iPCfe+KEdinXAmxpuUaVX7GmxS\n3NPNo0j90rrOGYMY65OHN+Fu2VuN78QRFiwQ09TqJsaoLiJxp3AmpjRv1Fa0LxeNSqxloyVoLFK/\nRXTmBjHSa+pW4c5PCOdgZ4bBtTgTD4mRdKU4uuoll+NtsVbNRFHty8Saty0pV6/EPR3qLFK/RXQW\nDmKMCAdjruiZ9+EH7ffQoofXnbIA16bSU4qvqVmMCYejcQ95RAfvpPz6nTGIUV8vVopDgkdTefcL\nb+ulNg1XwfWiN9c37+uEV/mGGGn9Qk/qd7HwONNM4oRuGewSe3Q2UrtFpfVb/9/vXuFo3CAOjMfF\nXq/ZPqrfOBb34iKxHu4zPcLRawa9focMGTJkyJD/JX8DEq9C4wMoYF0AAAAASUVORK5CYII=\n",
"prompt_number": 10,
"text": [
"\n",
"3\u22c5y\u2081 3\u22c5y\u2082 y\u2083\n",
"\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 - \u2500\u2500\n",
" 8 4 8 "
]
}
],
"prompt_number": 10
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"soln.pop(x)\n",
"\n",
"soln[x] = (x2 + x3) / 2\n",
"y_est2 = sym.simplify(y.subs(soln))\n",
"print '\\ny(x1 + x2) / 2 ->\\n'\n",
"y_est2"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"y(x1 + x2) / 2 ->\n",
"\n"
]
},
{
"latex": [
"$$- \\frac{1}{8} y_{1} + \\frac{3}{4} y_{2} + \\frac{3}{8} y_{3}$$"
],
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAIkAAAAeCAYAAAAYRz0yAAAABHNCSVQICAgIfAhkiAAABCdJREFU\neJzt20mIXEUYwPHfxCWaRcGLIXGZjBdBjRBcUZPBDQ/BYIyCxgWdIAHFGLeDeogHcXAhgqKgQdqF\nGMQtYmJAD+NyiLjFHRckXtSgB02MaFTi4es2L8309Hud6uVJ/6HpV13Vr//zdb16X1X10KdPE/ZJ\neK7pWINN+DXhedvNvOrjJNyIbfi2q0b5KJ33UqzELgx21aQ4P+PK6vFF+F10+F6nrN6l7CTHYmr1\neDH+Uo5gl9W7lJ0kyxrc3m2JFmir977tOnHJmIszsQOruuxShFJ698pIMkdrF8A1eB/T0uoUohX3\nXvDOTa90kop8HqdgK2ZXy0eLv2FxW6zyUdHcvaPek9px0hLxNz7D99XykEgAN3fNKB8d9U7VSZbg\nkerxKK5r0G4QD2I9Lq2ruxavJ/LJy3t4XPiuwAgW4Ju6doPK6U3vuTflYeyP5fi4ru4drE30ORVp\nb3ud8qYH3Tt5uzkdb2MnzsNXmbqpIlN/o4M+eSmrNyV0n4HJmIV/cEGm7hyReB2T6LMq0l2NnfSm\nB92zU63j8CgGcgpsxrKcbeHH6vPF+A0bMnVniGXmzwucD57A8eO8foTY09g5Tt2ImC7mpR3elNu9\nMLtyPLJsxEt1r43hhUx5Op4TAWuFinxXYxH3PN4ni0RyJV4TG3BFqWju3o6Yz8PlYh/uaZxdq0iR\nkwzkeGQ5El9mypNFcN+slpfiJlyYyG8iirg3854mhvNVopM8Jr6cWV32zuNOdJhJWI11eFl1L6gb\n6yTf4ZBMeRQH2J1ArRZB7jWaeQ/hVhxVLW/EgTitU4IT0MwdhsXoTYxE+9UqUu3dzBbz9D9xKF7B\nhw3a3iA6wkNie/sE8fuTjxK57A3340UxI6inmfcnONXu33QcXn3+uk2uqWP+aeZ4kbhQt6cUvjdz\nPCB2JfMwSSRXT45TtzdL/JUW3jtfLHUP52g7kXeNp0SnK0pFPvd2xHwubha3yinZN6Rgofh9A/Hl\n7mjQ7hl7LugsFMPg3Yk8amzDHwXaHyxmGl80qC/qPYIfRMCLkte9HTH/APfhXbwl8YbhcjFluwe3\nGX9qBz+JJWKYKYbmJQ3adnKz8BaRzI0ZfyQp4r1AdBLivj+YyLGelDGfcMMwVU6yVmTLi0WytqlB\nu2U4UfTWGbhELA13k/Pxqri3NyKv93yRH6yvtqsFf0s63f9IGfO2bxhOF8neFHHl3IlftH4F1TYL\nd4lANNosTMFMXJEpj8mXk4zHkEj06tcrDmpdryGpYw6XiaWHFXge59Yq8q6uTsQisej1QOa1u0Tm\n/GyC87eTq8TVVVt8ul7MEtaJ0aBXKV3M54iMPMuoWMApG1u0PpJ0ko7GPMX/3WwVSd/VYrNoWKwz\nNJqz9yKH4Q6cJXKK7fbcMe01/g8x79OnT58+ffr0KcS/hD5Noa4qiRsAAAAASUVORK5CYII=\n",
"prompt_number": 11,
"text": [
"\n",
" y\u2081 3\u22c5y\u2082 3\u22c5y\u2083\n",
"- \u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\n",
" 8 4 8 "
]
}
],
"prompt_number": 11
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 4"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"t, c, b, g, h = sym.symbols('t, c, b, g, h')\n",
"\n",
"P = sym.Function('P')\n",
"m = sym.symbols('m', real=True)\n",
"sol = sym.dsolve(sym.Derivative(P(t), t) - m * (c + b * P(t) - g - h*P(t)), P(t))\n",
"\n",
"sol"
],
"language": "python",
"metadata": {},
"outputs": [
{
"latex": [
"$$\\operatorname{P}\\left(t\\right) = \\frac{- c + g + e^{C_{1} + b m t - h m t}}{b - h}$$"
],
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAAkCAYAAAApQLALAAAABHNCSVQICAgIfAhkiAAABdxJREFU\neJzt3HuMXVUVwOGPaYUCaSivFguhRVJ8NBblJeVVpqa8UoiIZkRIJFgpAZTWUkhEQxNjoiYYiYBo\nKyG8yvBKUARFbWrig4IEA5EEQwyP8KhYEAotEKD8sfZ470zPvWfOnFt7Z2b/kpt7zt5nr73mzj5r\nr732OofMeOQj21uBTtGzvRXI1OIIfB19+BK+iROG0e64grJJuB0nVtThU1gyzGuPxwMV5bftY+II\nhGW6g0swAT9oKrsS141Q3lvp+48V283H3ypcf19F+W37yAN4dDIHp+OYIeVXYdMIZe6E/XG0sOL3\nYSEexHqciV9iGl7DPan/PjyLKfhvSR+n4EWchAPwRN0+8gAenZyIO7ClqaxH/NNbMRsL0vERYjDA\ne7hGuBU34Q9YhtuwL9ZgI76N3yYZ+6e2D+Fp3InJWFTQ7/O4Px1/HJdjprhR7q/YR6YmF3WJzD5c\nMaTsVOySjk8TVqwV5xSUXYpDsKsYxBNxd6o7NtXDtZiFPTFP+KYzlRvD/fDTdPwtHIkZdfvIi7hq\n7NUlMm8X/7tv4HM4WUzDmzBVDNAdKspcg0NxXvocjLWpbi7uTccbhTXfIGaALcIdeLdE/nFNMnZK\nes6s20fVP3Ks0SPu7vfF9DtJw0oUsSJ92nE0zsDD+AT2xvk1ZPZgKd7Bq8KyXlmiwwqxmHupRf1Z\nuKVExqhgvFvg68S09GO8gi/XlHcQVorpfbWwrk/WlPkzoeNPcBf2qCmPMTJ4GZuLuAu0D9Q/jH4c\nLlby+6Xy3+OvTddNE9a5eZY6RljpATbie03n38eNqZyIFqysIfNjOFtY8LPE1NscNitiKj6KXnET\njUtW4zHhe7wpgs+/we9S+ToR8oDpBv8DitgLu28TTUfOcrGSr8KKkvo3xI1BrMpfVj7LtZPZp1qM\nddzR6sc9U2M6vVbEBU8S4Y05+BFuFYuIMzSC4AMciBfEKpNwxpfiQ51SvAM8L27OZk6vKfMlMWjh\nC2Jh9X4Nef80+LfdQYSqxvva5X+0cyGOTd9rC+r6RezwKmFphnKq8NXWp/MtuEEM4h+OQM9twWp8\nEovF6n0XjVXwSLlYuDCP4ysdkPeo8HsvFj76zviVwfHfTAtuE9ZjSkFdDzanTxF3Kd6SvB47dkS7\n7cNlFa59UMOd6JTMTAVeEP5uEUcKK/Bei/r1+G5B+VJ8tr5qXcly4bMS4bPH5CjPNqfVDzwLH8af\nWtQvEdZ3flNZn1jorRMr4d50fmHTNX83PKs0GtlHhNHOxtfEvn8d/zczDFr5wAPpdkMH8GR8R6S3\nrTN4hdyfPotT/QJbuxgbNMJWRVwvtjOrsESxn/7/ZlnT8c3bTYtxRqsBPC99L9QYzBPS9Q+IhIxf\n4+2Ctr0i1lrkH2/W3gc+t0TfTGYQ7Szwy9rvTG3Abum7mePxixZtdsN/KuiXybSlyAeekT5lic1P\niU2MZmaL3aZWbaendplMRyiywAMuw9qStmtEBtHjTWW9ImPoL+l8ivCbn0vncw3eWh3KSny6pN+h\nLNO4YXJ8NGOVGAizS66biKuHlPUL/3eAK0Q0g9g9+nknFMxkBihyIXqFn/qPkrbviihFc9Rggsie\nJ8Jlm8UjJESS9a0jVTSTacd0kajzhLC+b4tpefkw2n5VI5nnYPxZ5KteqnGD7K2R/JMp53yRlTa3\nhoxpIjPuXx3RKJOpwGT8W8xodVikfYL+qCdvdXYn88UM2GqrfricIB6SHLPkAdydLBCuXJ94amTO\nCGT0iJj87iJNdJXiF5pkMh3nSXwmHS8U70qoyuF4RuMRpAtECuyYYiw+UjTamSFCjuvS+T5iEczw\nH5ci3IdVIo8YDsMjHdU0kyng8yKpaYB+W78DYjisxVHpeGeRGjBV9z3aVYvsA3cfr2u8PmmW8H/L\nHqMfyq5iI+qhdH6ysOib8MUO6Ng11A3TZDrP02IRN108V7fI1glTZRwmfN+B1zFNEhtO04R1L3sJ\nSSaTyWQymUwmk+lSPgBz7CJxXnrmlAAAAABJRU5ErkJggg==\n",
"prompt_number": 12,
"text": [
"\n",
" C\u2081 + b\u22c5m\u22c5t - h\u22c5m\u22c5t\n",
" -c + g + \u212f \n",
"P(t) = \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n",
" b - h "
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem 5"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"b, c, g, h, m= (-1.2, 2, 0.1, 1.1, 1.1)\n",
"\n",
"\n",
"def p_t(t, p0=1.):\n",
" \"\"\"\n",
" I used the sympy solution to write this by hand.\n",
" I could have done it with sym.solve and soln.subs,\n",
" but this was easier.\n",
" \"\"\" \n",
" # Constant of integration\n",
" cc = c - g - (h - b) * p0\n",
" \n",
" term1 = (c - g) / (h - b)\n",
" term2 = cc * np.exp(m * (h - b) * t) / (h - b)\n",
" return term1 - term2\n",
"\n",
"\n",
"times = np.linspace(0, 1, 100)\n",
"pbar = np.ones_like(times) * (c - g) / (h - b)\n",
"plt.plot(times, p_t(times, 0.5), 'r', label=r'P(0) = 0.5')\n",
"plt.plot(times, p_t(times, 1.3), 'g', label=r'P(0) = 1.3')\n",
"plt.plot(times, pbar, 'b', label='Equilibrium value')\n",
"plt.legend(loc=0)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "pyout",
"prompt_number": 13,
"text": [
"Legend"
]
},
{
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdYFHf+B/A3CCplwaX3Ir1IEYwt6mo0dk00GsvZYrlf\nvMRoqpfkEnLp7YzG82KM5aKmqicqtkiyamyoKBYUFUGaSO8d5vfH6CIKiuzCLMv79TzfZ2bZ3ZkP\n8yzvHb5TvnqCIAggIiKdoC91AUREpDkMdSIiHcJQJyLSIQx1IiIdwlAnItIhDHUiIh2idqh//PHH\nCAgIQI8ePTBt2jRUVlZqoi4iImoBtUI9OTkZa9asQWxsLM6fP4/a2lr89NNPmqqNiIgekYE6bzYz\nM4OhoSHKysrQqVMnlJWVwdHRUVO1ERHRI1Ir1C0sLPDKK6/AxcUFRkZGGD58OIYOHap6Xk9PT+0C\niYg6opZe7K9W90tiYiK++uorJCcnIyMjAyUlJdi8efN9hbEJePfddyWvQVsatwW3BbfFg5s61Ar1\nU6dOoV+/frC0tISBgQEmTJiAo0ePqlUQERG1nFqh7uvri+PHj6O8vByCIODAgQPw9/fXVG1ERPSI\n1Ar14OBgzJw5E+Hh4QgKCgIALFiwQCOF6RqFQiF1CVqD26Iet0U9bgvN0BPU7cB50ML19NTuHyIi\n6mjUyU61zn5pKQsLC+Tn50uxatJBcrkceXl5UpdBpBUk2VPnHjxpEj9PpGvU+Uzz3i9ERDqEoU5E\npEMY6kREOoSh3gKrV6/GkiVLmvXalStXYunSpa1cERGRiKHeCDc3NxgbG0Mmk8HOzg5z5sxBaWkp\nAKCqqgoffvghXn/9ddXrz549i7CwMJiYmCA8PBxxcXGq5+bPn4/NmzcjOztb43U+aL33mj17Nrp0\n6QKZTAaZTAYzMzMeXCTSQQz1Rujp6WHXrl0oLi5GbGwsTp06hQ8++AAAEBkZCT8/P9jb2wMQQ378\n+PGYOXMmCgoKMGvWLIwfPx7V1dUAgC5dumDkyJH4/vvvNVrjw9bb2O/0xhtvoLi4GMXFxSgqKuIN\n14h0EEP9IRwcHDBixAhcvHgRALBnzx4MGjRI9bxSqURtbS1eeuklGBoa4sUXX4QgCPj9999Vr1Eo\nFIiKitJoXc1Z7724Z06k+xjqTbgTgKmpqdizZw9CQ0MBABcuXICPj4/qdRcvXlTdIuGO4OBg1ZcA\nIN4j50FdI0FBQZDL5Y22F154odH3NGe991q1ahUsLS0RHh6Obdu2Nfk6Imq/JLmitFk00TXQwj1T\nQRDw1FNPwcDAAObm5hgzZgzefPNNAEBBQQFkMpnqtSUlJTA3N2/wfjMzMxQXF6sey2QyFBYWNrm+\nc+fOPXKNzVnv3RYtWoR//etfMDc3x759+/Dss8/Czs4O/fr1e+R1E5H20t5Ql7CrQE9PD5GRkRgy\nZMh9z8nlchQVFakey2SyBo8BoLCwEGZmZqrHxcXF9wWwupqz3rvd+U8DAEaOHInp06dj27ZtDHUi\nHcPul0cUFBSEK1euqB4HBATct6d97tw5BAQEqB5funQJISEhTS4zICBAdVbKvW3hwoVNvudh6yWi\njoeh/ohGjRqFgwcPqh4rFAp06tQJK1asQGVlJVasWAF9ff0Ge/kHDx7EyJEjm1zmxYsXVWel3NtW\nrVrV6Huas967bdmyBSUlJairq8P+/fuxefNmjBs3roVbgYi0ltCKmlp8K69WbW5ubkJ0dHSjz1VV\nVQkuLi5CRkaG6mdnzpwRwsLCBCMjIyEsLEw4e/as6rny8nLByclJyMrK0nidD1rvpk2bhICAANXj\nAQMGCObm5oKZmZkQEhIi/PzzzxqvRyra/nkielTqfKZ5l8YWWLNmDeLj47Fs2bKHvnblypVIS0vD\nJ5980gaVdUzt/fNEdC91PtMMdWr3+HkiXcNb7xIR6YgbBTfUej9DnYhISxy+cRh91vZRaxkMdSIi\nLbA2di0m/jIRG8ZvUGs52nvxERFRB1BdW40l+5bgt+u/4fCcw/Cx8nn4mx6AoU5EJJGs0ixM+nUS\nZJ1liJkXA/Ou6l95zu4XIiIJnMo4hV5remGAywDsmLpDI4EOcE+diKjNrY1di6XRS/HN6G8w0X+i\nRpfNUG+B1atX4/Lly7z4iIgeSUVNBRbtWYTDKYdxaPYh+Fn7aXwd7H5pRHsZzm7BggXw9fVFp06d\n8N///veBr3399dfh4uICMzMzODk54eWXX0ZNTY3GayKixiUXJGPA+gHIK89DzLyYVgl0gKHeqPYw\nnB0AhISEYNWqVejZs+dDh6abO3cu4uPjUVRUhJiYGOzfvx/fffedxmsiovvturILvb/rjamBU/Hr\npF8h6yJ7+JtaiKH+ENo6nB0ALFy4EEOGDEHXrl0f+lofHx+YmpoCEAcB0dfXV30xEVHrqKmrwZvR\nb+L5qOexbfI2vNz35VYfG1jtUC8oKMAzzzwDPz8/+Pv74/jx45qoS3KClg9n1xKffPIJZDIZnJ2d\nMWbMGIwfP15jyyaihtKK0jDkv0NwKuMUTi84jf4u/dtkvWofKH3ppZcwatQobNmyBTU1Naq+Z3Xp\nvaf+t5nwru4OZ9cSS5cuxdKlS3HmzBk89dRTCA8Px4QJE9pk3UQdye6ru/Fc5HNY1HsRlj6+FPp6\nbdcpolaoFxYW4vDhw6qDdHdCUBNaGsia0B6Gs1NHaGgoFi5ciI0bNzLUiTSoqrYKb//+Nn688CN+\nnfQrBrgOaPMa1Ar1pKQkWFtbY86cOYiLi0NYWBiWL18OY2Nj1WsiIiJU8wqFAgqFQp1VSq6x4ey+\n/PLLBq85d+4cXnzxRdXj5gxnl5KS0uhzM2bMaHL0I3VUV1fDxMRE48sl6qiu5V3D1K1TYW9qjzN/\nPQMrY6tmv1epVEKpVGqmEDUG5xBOnjwpGBgYCDExMYIgCMJLL70k/OMf/1A939Ti1Vxtq3vQyEfb\ntm0TnnzySdXjqqoqwdXVVVi+fLlQUVEhLF++XHBzcxOqq6tVr5k/f77w+eefa7zOqqoqoby8XOjX\nr5+wZs0aoby8XKirq7vvdXV1dcI333wj5OfnC3V1dcKJEycEe3t7YevWrRqvSQra/nki3VZXVyd8\nf/Z7weozK+HrE183+jf4qNT5TKv113Dz5k3Bzc1N9fjw4cPC6NGjH1qYtv8Rtpfh7AYNGiTo6ekJ\n+vr6gp6enqCnpyccPHhQEISGw9nV1tYKI0aMECwsLASZTCYEBgYKa9eu1Xg9UtH2zxPprvzyfGHK\nlimC30o/4ezNsw9/QzOp85lWe+SjgQMH4rvvvoO3tzciIiJQXl6OTz/9FIDujnzE4ey0S3v/PFH7\ndDD5IGZun4lxPuPw2dDPYGRopLFlSzqcXVxcHObNm4eqqip4eHhg/fr1qoOCuhrqpF34eaK2VFlT\niX/88Q9sOrcJ3437DqO8Rml8HRyjlDo0fp6orcRlxmHG/2bA08ITq8eshrWJdausR53PNG/oRUT0\nEDV1Nfji6Bf417F/4Ysnv8CMoBmtfmVoSzHUiYge4HLOZczePhumnU1xcv5JuHZzlbqkB+K9X4iI\nGlFbV4svj36JAesHYFbwLPw24zetD3SAe+pERPeJz47Hc5HPoatBV5yYdwLd5d2lLqnZuKdORHRb\ndW01Pjr8EQauH4hZwbPw+6zf21WgA9xTJyICII4ZOm/HPNia2uL0gtPtoqulMdxTb0MpKSmQyWSq\nU5UUCgXWrl0LANi8eTOGDx+ueq2+vj6uX7/e7GXf+35tlJycDH19fdTV1UldCpFKaVUpXtn/Csb8\nMAav9nsVe6fvbbeBDjDUG3X3cHZ32qJFi9RerouLC4qLi1WnQunp6anmp0+fjn379rV42eq+n6gj\n2n11NwL/E4hbJbdw/vnz+EvQX7T2VMXmYvdLI+4MZ9fYrXe1UW1tLTp16iR1GUTtRkZxBhbvXYzY\nm7FYPWY1nvR4UuqSNIZ76o+orq4Or776KqytreHh4YF///vfDboU3NzcEB0drXp9REQEZsyYAeDB\n3Q8bNmzAgAEN770cFRUFDw8PWFtb4/XXX1d122zYsAH9+/fHyy+/DCsrK0RERDR4f2Prubur5+73\ny+VyeHp64ujRo1i/fj1cXFxga2vb5JiqP//8M3r16tXgZ8uWLVONohQVFYXQ0FCYm5vDxcUF7733\nXpPb8kHbCgCOHz+Ofv36QS6XIyQkBAcPHmxyWUTNUVNXgxUnViD4m2B4W3rj/PPndSrQAYZ6k5q6\nRPfbb79FVFQUzp49i1OnTmHLli0N/l27u0vlzuOW2r59O06fPo3Y2FhERkZi3bp1qudiYmLg4eGB\nrKwsvPXWWw9d1r11xcTEIDg4GHl5eZg6dSomT56M2NhYJCYmYtOmTXjhhRdQVlZ233LGjRuHhIQE\nXLt2TfWzH374AdOnTwcAmJqaYtOmTSgsLERUVBT+85//IDIyslk13T2fnp6OMWPG4J133kF+fj6+\n+OILTJw4ETk5OQ/9XYkaczztOHqt6YXtl7fj0OxD+GDIBxq9CZe20NpQ19NTv7WUcHs4u7vHCr2z\nl/vLL79gyZIlcHR0hFwux5tvvvnAezSoc0+SN954A926dYOzszMWL16MH3/8UfWcg4MD/va3v0Ff\nX79ZA0/fy93dHbNmzYKenh4mT56MjIwMvPPOOzA0NMSwYcPQuXPnBsF9h5GREcaPH6+q5erVq0hI\nSMC4ceMAAIMGDUJAQAAAoEePHpgyZUqz97Dv3labNm3CqFGjMGLECADA0KFDER4ejt27dz/y70od\nW3ZpNubtmIeJv0zEa/1eQ/TMaPhZ+0ldVqvR2lAXBPVbS90Zzi4/P1/V5s6dCwC4efMmnJ2dVa91\ncXFR91dt0r3rycjIaPS5lrC1tVXNGxmJeyvW1tYNflZSUtLoe6dNm6YK9R9++AFPP/206ovlxIkT\nGDx4MGxsbNCtWzesXr0aubm5j1zfjRs38Ouvvzb4Yj1y5AgyMzMfeVnUMdXU1WBlzEr4r/KHWRcz\nxC+Mx7Qe09r9gdCH4YHSR2Rvb99g6Ll7h6EzMTFpMPi2OiGUkpICPz8/1byjo6PquQd9MO8MU1dW\nVgZTU1O167jX0KFDkZ2djbi4OPz000/46quvVM9NmzYNixYtwr59+9C5c2csWbKkyS6TxrbVnd/L\nxcUFM2bMwLfffquxuqnj+CPpDyzetxhWxlZQzlIiwCZA6pLajNbuqUutqW6TyZMnY8WKFUhPT0d+\nfj4++eSTBgEbEhKCn376CTU1NTh16hS2bt3a4j2DL774AgUFBUhNTcWKFSvw7LPPNut91tbWcHR0\nxMaNG1FbW4t169YhMTGxRTU0xtDQEJMmTcKrr76K/Px8DBs2TPVcSUkJ5HI5OnfujJiYGPzwww9N\n/v6Nbas7/vKXv2Dnzp3Yv38/amtrUVFRAaVSifT0dI39HqR7rudfx8RfJuK5Hc/hHwP/gQMzDnSo\nQAcY6k0aO3Zsg/PUJ06cCACYP38+hg8fjuDgYISHh2PixIkNvgDef/99JCYmQi6XIyIiQnUA8Y6m\nAu7eg4YAMH78eISFhSE0NBRjxoxRdQE19tp7f7ZmzRp8/vnnsLKyQnx8PPr37//AdT3qF8+0adMQ\nHR2NSZMmQV+//mO0atUqvPPOOzAzM8P7779/3xfR3et50LZycnJCZGQkPvroI9jY2MDFxQVffvkl\nL1yiRhVVFuHv0X/HY2seQ0+7nohfGI9n/J/R+a6WxnCQDDUlJyeje/fuqKmpaRBu1HZ06fNEj6am\nrgZrY9ci4mAERnqOxAdDPoCDzEHqstTGQTKIqEMRBAFRV6PwxoE3YGNig93TdiPUPlTqsrQCQ10D\nOuK/eERSOZVxCq/99hpuldzCZ8M+w2iv0fwbvAu7X6jd4+epY7iaexVv/f4W/kz5ExGKCDwX+hwM\n9HVzv1SdzzQ7gYlIq90svonno55H37V9EWIXgqsvXsWCsAU6G+jq4lYhIq2UW5aLT498irVn1mJO\nyBwkvJAAS2NLqcvSegx1ItIqhRWFWHZ8GVbGrMTkgMk493/n4Gjm+PA3EgCJQl0ul/PABmmMXC6X\nugTSgOLKYqw4sQJfnfgKo7xGIWZ+TLsbSk4bSBLqeXl5UqyWiLRQcWUxVsasxFcnvsKw7sNw5Lkj\n8Lb0lrqsdovdL0QkicKKQnwd8zVWnFiBYR7DoJyl1Om7J7YVhjoRtancslwsP7Ecq06uwkivkTg8\n5zB8rHykLktnMNSJqE1kFGdg2fFlWHdmHSb4TcDxecfhaeEpdVk6R+3z1GtraxEaGoqxY8dqoh4i\n0jHX8q5hwc4FCFwViKraKpz961msGbuGgd5K1N5TX758Ofz9/VFcXKyJeohIR8Skx+Dzo59DmazE\n8+HP48qLV2BlbCV1WTpPrT31tLQ07N69G/PmzeNl2kSEOqEOUVeioNigwKRfJ6G/c39cX3Qd/xz8\nTwZ6G1FrT33JkiX4/PPPUVRU1ORrIiIiVPMKhQIKhUKdVRKRFiqvLsfGcxux7PgydDXoilf6voJn\nA56FYSdDqUtrF5RKJZRKpUaW1eIbeu3atQt79uzBv//9byiVSnz55ZfYuXNnw4XzRktEOi2jOAP/\nOfUffHv6Wzzm+Bhe7vMyFG4KXlyoJknup3706FHs2LEDu3fvRkVFBYqKijBz5kx8//33LV0kEbUT\nJ9NPYkXMCkRdicLUHlNxaPYhnpaoJTRy692DBw/iiy++4J46kQ6rrKnEr/G/4uuYr5FVmoW/9fob\n5obOhdyIt2nQNK0Y+Yj/bhHpppTCFKw+vRprY9eih20PvDXgLYz2Go1O+p2kLo0aIckgGUSk3Wrr\narE/cT++Of0N/kz5EzOCZuD/wv8Pvla+UpfWIaiTnQx1IlLJLMnEujPr8O3pb2FlbIW/hv0V03pM\ng0lnE6lL61C0ovuFiNqn2rpa7EvchzWxa6BMVuIZ/2ewZfIWhDuES10atQD31Ik6qOv517Hh7AZs\nOLsBdqZ2mN9zPqYEToGsi0zq0jo87qkTUbOUVpVi26VtWH92Pc5nncf0HtOxc+pOBNsFS10aaQj3\n1Il0XJ1Qh8M3DuO/cf/F/y7/D/2c+2F28GyM8xmHLgZdpC6PGsEDpUR0n4ScBGw8txGbzm2CaWdT\nzA6Zjek9psNeZi91afQQ7H4hIgDi2Ss/XfgJm89vRmphKqb1mIbtU7Yj2DaY15J0ENxTJ2rn8svz\n8b/L/8OPF37EqYxTGO8zHtN6TMMQ9yEw0Od+W3vE7heiDqa4shg7Enbgl/hfoExWYmj3oZgSMAWj\nvUfD2NBY6vJITQx1og6guLIYUVej8MvFXxCdFI0BLgMwyX8SnvJ9CuZdzaUujzSIoU6kowoqCrDr\nyi5sid+CP5L/wOMuj2Oi30Q87fs0b6SlwxjqRDoksyQTOxJ2YNulbTiaehSD3Qdjgu8EjPMZxyDv\nIBjqRO1cQk4CIhMisf3ydlzKuYQRniMwwXcCRnqNhGlnU6nLozbGUCdqZ2rqanAs9Rh2XtmJyIRI\nlFSVYKz3WDzt+zQUbgpeFNTBMdSJ2oH88nzsS9yHXVd2Yc+1PXA2c8Z43/EY5z0OPe178jxyUmGo\nE2khQRAQdysOe67uwe5ruxGXGYeBrgMx1nssRnuPhpOZk9QlkpZiqBNpidyyXBy4fgB7E/di77W9\nMO1sihGeIzDaazQGuQ6CkaGR1CVSO8BQJ5JIVW0VjqUew2/Xf8P+xP1IyE3AQNeBGOExAsM9h8PT\nwlPqEqkdYqgTtZE6oQ4Xsi4g+no0DiQdwOEbh+Ft6Y1hHsMw3GM4+jn3Q+dOnaUuk9o5hjpRKxEE\nAdfyruGP5D/we9Lv+D3pd5h3NccQ9yEY6j4UQ9yHwNLYUuoySccw1Ik0RBAEJOYnQpmsxMEbB/FH\n0h/Q09PDYLfBGOw2GEPch8C1m6vUZZKOY6gTtVCdUIdL2ZdwOOUwDt04hEM3DgEAFG4KDHIdBIWb\nAp4WnjzdkNoUQ52omSprKhF7MxZ/pvyJwymHcST1CLp17YaBrgMx0GUgBroORHd5d4Y4SYqhTtSE\nrNIsHEs9hqNpR3Ek5QjOZJ6Bt6U3Hnd5HANcBuBxl8fhIHOQukyiBhjqRACqa6tx7tY5HE87juPp\nx3Es9Rhyy3PR27E3+jr1RX+X/ujt2BuyLjKpSyV6IIY6dTiCIOB6/nWczDiJmPQYnEg/gbjMOLjL\n3dHHqQ/6OPZBH6c+8LP2g76evtTlEj0ShjrpNEEQkF6cjtMZp3Hq5imcTD+JUxmnYGRohF4OvdDb\nsTcec3wMYQ5hMOtiJnW5RGpjqJPOEAQBqUWpiL0Zi9M3T4vTjNOoE+oQ7hCuar0cesFeZi91uUSt\nQrJQT01NxcyZM5GVlQU9PT0sWLAAixYt0khhpPtq6mqQkJOAs5lncfbWWZy5eQZnMs+gc6fOCLUL\nRZhDGMLsw9DTvieczZx5Rgp1GJKFemZmJjIzMxESEoKSkhKEhYVh+/bt8PPzU7sw0i05ZTk4f+s8\nzt06h7hbcYi7FYdL2ZfgbO6MELsQBNsGI9QuFKH2obAztZO6XCJJqZOdBuqs2M7ODnZ24h+gqakp\n/Pz8kJGRoQp16njKq8txKecSzt86j/NZt9ut8yitLkUPmx4Isg1Cb8feWBC2AIE2gRzVh0jD1Ar1\nuyUnJ+PMmTPo3bt3g59HRESo5hUKBRQKhaZWSRIqry7H5ZzLiM+OR3xOPC5mXcSFrAtIL06Hl4UX\netj2QKB1IF7o9QKCbIPgYu7C7hOiJiiVSiiVSo0sSyMHSktKSqBQKPD222/jqaeeql84u1/avbzy\nPFzOuYzLOZdxKecSLmVfQnx2PDKKM+Bl6QV/a3/4W/nD39ofgTaB8LTwhGEnQ6nLJmrXJD37pbq6\nGmPGjMHIkSOxePFijRVGbaeqtgrX86/jSu4VJOQkICH3dstJQEVNBfys/eBr5QtfS1/4WfvBz8oP\n3eXdGd5ErUSyUBcEAbNmzYKlpSWWLVum0cJIs6prq3Gj8Aau5l7FtbxruJp3FVfzruJK7hWkF6XD\nxdwF3pbe8LL0go+lj9isfGBvas9uE6I2Jlmo//nnnxg4cCCCgoJUf/gff/wxRowYoXZh9OhKqkpw\nPf86rudfR2JeIhLzxXYt7xrSitLgIHOAl4UXPC084WXhBS9LL3hZeMFd7s6BHYi0CC8+6iAqayqR\nWpSKpPwkJBUkIbkgGUkFSbiefx1J+UkoqSqBu9wdHnIPdJd3h4fcAx4WHvCQe8Ctmxu6GHSR+lcg\nomZgqOuIkqoSpBSmIKUwBTcKbuBG4e1WcAPJBcnILsuGg8wB7t3c4S53h5u5G9zl7ugu7w73bu6w\nM7VjVwmRDmCotwOVNZXIKM5AalEqUgtTxent+ZTCFKQWpaKsugwu5i5wNXdVTV27ucLV3BXucnc4\nyBxgoK+xs1CJSEsx1CUkCAIKKwuRXpSOjOIMpBenI70oXZwWpyOtKA3pRenIK8+DvcwezmbOcDJz\ngrO5M5zNnOFi7gJnM2c4mzvD2tiae9pExFBvDYIgIK88D5klmbhZchM3i2+K09vzGcUZqtZJvxMc\nZY5wNHOEo8wRDjIHOJk5wcnMCY4yRziZOcHGxAad9DtJ/WsRUTvAUG8mQRCQX5GPrNIs3Cq5hVul\nt+qnt+czSzKRWZKJW6W3YGRgBHuZPexN7WFnagcHmQPsTe1VP3M0c4S9qT0HXSAijeqwoX6n6yOn\nLAfZpdnILstGdmk2skqzxPmy2/Ol2bhVegvZpdkwNjSGjYkNbE1tYWtiCxsTG9iZ2sHWxFacmtrC\n3tQetqa26GrQtdVqJyJqik6Eep1Qh4KKAuSW5SKvPA+55bnIKctBbllu/fzt6Z2WW5aLrgZdYW1i\nDWtja1gZW8HGxEb12MbERnx81zxP6yMibafVoR55ORJ55XnIL89HfkW+OH97encrrCiErIsMlkaW\nsDCygKWxJSyNLGFlbFU/vf0zaxNr1c8Z0kSka7Q61EdvHg0LIwvIjeSQd5XDwshCfHx73tJYDPFu\nXbvxdD0iImh5qGvTgVIiovZAnezkMOtERDqEoU5EpEMY6kREOoShTkSkQxjqREQ6hKFORKRDWv3E\ncN50kIio7bR6qPM0dSKiR6POzjC7X4iIdAhDnYhIhzDUiYh0CEOdiEiHtP5tEdeuBUxMAFPThk0m\nE6cmJoAB785IRKQJrZ+mR48CJSViKy4GSkvr5+/8vEsXMeQf1MzM6qcPavyCIKIOTPpb7woCUFZW\nH/LFxfe3oqL66d3zhYUNf15UBHTtKoa7ufn9rVu3hvONNZmMJ9cTkaR4P/U7BEH8T6CwUGx3gv/e\nlp8vTgsKxHb3fFlZfejL5WK7M29hcf/0TrO0BIyN+YVARGpjqGtSTU198DfW8vLun+bmivM1NQ1D\n/u5mZVU/vXteLgf0ebyaiOox1LVFRUV9yDfWcnIazmdni11OcrkY8NbWYrOyAmxsxPm7pzY24pdB\np05S/6ZE1IokDfW9e/di8eLFqK2txbx58/DGG29opLAOo7pa/CLIzq4P+rtbVpbY7swXFIhfAra2\nYsjb2jZsdnb1UxsbHjgmaockC/Xa2lr4+PjgwIEDcHR0RK9evfDjjz/Cz89P7cKoCTU1YvhnZQG3\nbjXeMjPFlpsrfgHY24shb2/fsDk41E+7dpX6NyOi29TJTrV242JiYuDp6Qk3NzcAwJQpUxAZGakK\ndWoFBgZiQNvZPfy1tbXiHn5mJnDzZn1LSACUSnE+I0N83tRUDHcHB8DRsWFzchKblRUPBBNpObVC\nPT09Hc7OzqrHTk5OOHHiRIPXREREqOYVCgUUCoU6q6RH0alT/RdASEjTr6urE/fqMzKA9PT6aWws\nsHMnkJYmPi4pEUPe2bk+6J2dxebiIja5nMFP9IiUSiWUSqVGlqVWqOs144/37lAnLaWvX3+QNji4\n6deVl4sBn5YGpKaK7fx5YPducT4lReweuhPwrq71U1dXwM1N/E+AB3qJGrh3h/e9995r8bLUCnVH\nR0ekpqYH6MTAAAANcUlEQVSqHqempsLJyUmdRZI2MzICvLzE1pSiIjHcU1KAGzfEFhUlTpOTxf8I\nnJzEgHdzA9zd61v37uJBXu7pE7WYWgdKa2pq4OPjg+joaDg4OOCxxx7jgVJ6sIqK+sBPSqpv16+L\n07Ky+oD38KifeniIP+/cWerfgKjVSXpK4549e1SnNM6dOxd///vfNVIYdVBFRWK4JyaKQX/9OnDt\nmvg4LU3svvH0FP9buDP19mbgk07hxUfUMVRXi104iYnA1ati2F+9KrbUVLFbx9tbbD4+4tTXV/wi\nYJcOtSMMdaKqKnEP/8oVsSUk1LfSUjHkfX0BPz9x6u8v7ukbGkpdOdF9GOpED5KfD1y+LAb8pUv1\nLTVV7LP39xdbQIA49fFhVw5JiqFO1BIVFWLXzcWLQHy82C5eFLt4uncXQz4wEOjRQ2zdu/Pma9Qm\nGOpEmlRZKe7Vnz8PXLggTs+fF0/HDAgQz+UPChJbcLB4q2YiDWKoE7WFwkIx3M+dA+LixOn58/UX\nbYWGilfuhoaKV9ny4Cy1EEOdSCp1deLZOGfPiu3MGbFVVQE9e9a3sDDxXHsGPTUDQ51I29y8KYZ7\nbKzYTp8W9/R79gR69QLCw8Xm5sagp/sw1Inag+xsMdxPnwZOnhRbVZUY8r16AY89BvTuLd4Nkzo0\nhjpRe5WeXh/wJ06IU2trMdz79AH69hX763k+fYfCUCfSFbW14pk3x4/Xt+vXxW6bvn2Bfv3EZm0t\ndaXUihjqRLqsqAiIiQGOHhXb8ePi3Sz79wcef1xsXl7sm9chDHWijqS2VrxI6sgR4M8/gcOHxb75\nxx8HBg4UW1AQL5RqxxjqRB3djRvAoUNiwB88KI5h+/jjgEIBDB4s9stzcJJ2g6FORA1lZoohr1TW\nj0c7aBAwZIjYAgLYXaPFGOpE9GCZmWK4//478Mcf4jnzTzwhtqFDxfPlSWsw1Ino0dy4AURHi+3A\nAcDMDBg2TGxDhvB+NhJjqBNRy9XVifew+e03YP9+4Ngx8R42w4eLLSyMB13bGEOdiDSnvFzsj9+3\nD9i7F8jJEcN95EhxamkpdYU6j6FORK3nxg0x3HfvFvvlAwOB0aPFFhTEA66tgKFORG2jslI8ZTIq\nCti1C6ipAcaMAcaOFU+d7NJF6gp1AkOdiNqeIIjDAu7cKbYLF8QDrePGiXvxFhZSV9huMdSJSHpZ\nWeIefGSkeOrkY48BTz8NPPUU4OgodXXtCkOdiLRLaal4Js3//id20/j4ABMmABMnimO90gMx1IlI\ne1VViRc8bd0KbN8uDvU3aZLYPDykrk4rMdSJqH2oqRFPl/z1V2DbNsDJCXj2WWDyZF7VeheGOhG1\nP7W14pk0P/8s7sV7eQFTp4oBb2cndXWSYqgTUftWXS3eruDHH4EdO8Th/aZPF/vhzcykrq7NMdSJ\nSHeUlYln0WzeLPbFjxgBzJghXs3aQYb1Y6gTkW7KzRX73zduBK5dE7tnZs4EQkN1+kpWSUL9tdde\nw65du9C5c2d4eHhg/fr1ML/nzm4MdSLSmGvXxHD//ntAJgPmzBG7aGxspK5M49TJzhbfeu3JJ5/E\nxYsXERcXB29vb3z88cctXRQR0cN5egLvvQckJgLLlwNnzgDe3uIFTnduWUAtD/Vhw4ZB//btOHv3\n7o20tDSNFUVE1CR9ffE+M99/D6SkAKNGAR98ALi6Am+9BSQlSV2hpAw0sZB169Zh6tSpjT4XERGh\nmlcoFFAoFJpYJRGReGbM/Pliu3gR+O478fYEoaHAggXifWg6d5a6yodSKpVQKpUaWdYD+9SHDRuG\nzMzM+37+0UcfYezYsQCADz/8ELGxsdi6dev9C2efOhG1tYoK8cKm1auBK1eAuXPFgHdxkbqyZpPs\n7JcNGzZgzZo1iI6ORteuXTVaGBGR2i5dAr75Bti0CejXD/jb34Ann9T6kZwkCfW9e/filVdewcGD\nB2FlZaXxwoiINKasTLywaeVKoKQEeP554LnngG7dpK6sUZKEupeXF6qqqmBx+57Jffv2xapVqzRW\nGBGRxgmCOAbrypXAnj3iee8vvAD4+0tdWQO8+IiI6FFlZIj97qtXA8HBwOLF4lWrWtA1w1AnImqp\nykrgp5+AZcvE+ZdeAmbNAoyMJCuJoU5EpC5BEAfWXrYMOHFC7HdfuFCSK1YluaKUiEin6OmJFzXt\n2CHeEvjmTXHEpv/7P/EWBe0EQ52I6F6+vmJf++XLgJUV0KePeJ/32FipK3sohjoRUVNsbcVbECQl\nicE+bpx4K+BDh8TuGi3EPnUiouaqrBQvZPrkEzHw335bPGNGw7cB5oFSIqK2VFsr3uf9gw/Es2Te\nflvci9dQuDPUiYikUFcHbN8OvP+++Pjdd4Hx49UOd4Y6EZGUBAHYuROIiBCDPiJCrXBnqBMRaQNB\nEE+JfPddwMAA+Oc/gZEjHzncGepERNqkrk68/e+77wLm5sCHH4rnwDcTQ52ISBvV1oq3IHjnHcDD\nA/joIyA8/KFv4xWlRETaqFMncXDsy5eBCRPEfvZJk8TBO1oJQ52IqLUZGoq3G7h6FejZUxywY+FC\noJGR5dTFUCciaivGxsDf/w4kJABduwIBAeLB1NJSja2CoU5E1NYsLYF//Qs4dUoccs/bG1i7VuyD\nVxMPlBIRSS0mBnjlFaC4GFi2DHpDhvDsFyKidk0QgK1bgddeg15yMkOdiEgnVFRAz8iIoU5EpCt4\nnjoREQFgqBMR6RSGOhGRDmGoExHpEIY6EZEOYagTEekQhjoRkQ5hqLcRpVIpdQlag9uiHrdFPW4L\nzVA71L/88kvo6+sjLy9PE/XoLH5g63Fb1OO2qMdtoRlqhXpqaip+++03uLq6aqoeIiJSg1qh/vLL\nL+Ozzz7TVC1ERKSmFt/7JTIyEkqlEsuWLYO7uztOnz4NCwuLhgt/xBG0iYhI1NJ7vxg86Mlhw4Yh\ns5Hhlj788EN8/PHH2L9//wML4M28iIjaVov21C9cuIAnnngCxsbGAIC0tDQ4OjoiJiYGNjY2Gi+S\niIiaRyO33m2q+4WIiNqWRs5TZ985EZF20Eior1q1Cv369YOXlxc+/fTTRl+zaNEieHl5ITg4GGfO\nnNHEarXS3r174evr2+S22Lx5M4KDgxEUFIT+/fvj3LlzElTZNh62Le44efIkDAwMsG3btjasrm01\nZ1solUqEhoYiMDAQCoWibQtsQw/bFjk5ORgxYgRCQkIQGBiIDRs2tH2RbeC5556Dra0tevTo0eRr\nWpSbgppqamoEDw8PISkpSaiqqhKCg4OF+Pj4Bq+JiooSRo4cKQiCIBw/flzo3bu3uqvVSs3ZFkeP\nHhUKCgoEQRCEPXv2dOhtced1gwcPFkaPHi1s2bJFgkpbX3O2RX5+vuDv7y+kpqYKgiAI2dnZUpTa\n6pqzLd59911h6dKlgiCI28HCwkKorq6WotxWdejQISE2NlYIDAxs9PmW5qbae+oxMTHw9PSEm5sb\nDA0NMWXKFERGRjZ4zY4dOzBr1iwAQO/evVFQUIBbt26pu2qt05xt0bdvX5ibmwMQt0VaWpoUpba6\n5mwLAPj666/xzDPPwNraWoIq20ZztsUPP/yAiRMnwsnJCQBgZWUlRamtrjnbwt7eHkVFRQCAoqIi\nWFpawsDggSfqtUsDBgyAXC5v8vmW5qbaoZ6eng5nZ2fVYycnJ6Snpz/0NboYZs3ZFndbu3YtRo0a\n1Raltbnmfi4iIyPx/PPPA9DdYzPN2RZXr15FXl4eBg8ejPDwcGzcuLGty2wTzdkW8+fPx8WLF+Hg\n4IDg4GAsX768rcvUCi3NTbW//pr7hyjcc5KNLv4BP8rv9Mcff2DdunU4cuRIK1YkneZsi8WLF+OT\nTz5RDbJ772dEVzRnW1RXVyM2NhbR0dEoKytD37590adPH3h5ebVBhW2nOdvio48+QkhICJRKJRIT\nEzFs2DDExcVBJpO1QYXapSW5qXaoOzo6IjU1VfU4NTVV9S9kU6+5c167rmnOtgCAc+fOYf78+di7\nd+8D//1qz5qzLU6fPo0pU6YAEA+O7dmzB4aGhhg3blyb1tramrMtnJ2dYWVlBSMjIxgZGWHgwIGI\ni4vTuVBvzrY4evQo3nrrLQCAh4cH3N3dkZCQgPDw8DatVWotzk11O/urq6uF7t27C0lJSUJlZeVD\nD5QeO3ZMZw8ONmdb3LhxQ/Dw8BCOHTsmUZVtoznb4m6zZ88Wtm7d2oYVtp3mbItLly4JTzzxhFBT\nUyOUlpYKgYGBwsWLFyWquPU0Z1ssWbJEiIiIEARBEDIzMwVHR0chNzdXinJbXVJSUrMOlD5Kbqq9\np25gYICVK1di+PDhqK2txdy5c+Hn54fVq1cDAP76179i1KhR2L17Nzw9PWFiYoL169eru1qt1Jxt\n8c9//hP5+fmqfmRDQ0PExMRIWXaraM626Ciasy18fX0xYsQIBAUFQV9fH/Pnz4e/v7/ElWtec7bF\nm2++iTlz5iA4OBh1dXX47LPPdPLCxqlTp+LgwYPIycmBs7Mz3nvvPVRXVwNQLzc1ckUpERFpB458\nRESkQxjqREQ6hKFORKRDGOpERDqEoU5EpEMY6kREOuT/AXJLVglb8EqlAAAAAElFTkSuQmCC\n"
}
],
"prompt_number": 13
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}