eg2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "import numpy as np\n",
    "from matplotlib import pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def bisection(f,a,b,N):\n",
    "    '''\n",
    "    https://www.math.ubc.ca/~pwalls/math-python/roots-optimization/bisection/\n",
    "    Approximate solution of f(x)=0 on interval [a,b] by the bisection method.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    f : function\n",
    "        The function for which we are trying to approximate a solution f(x)=0.\n",
    "    a,b : numbers\n",
    "        The interval in which to search for a solution. The function returns\n",
    "        None if f(a)*f(b) >= 0 since a solution is not guaranteed.\n",
    "    N : (positive) integer\n",
    "        The number of iterations to implement.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    x_N : number\n",
    "        The midpoint of the Nth interval computed by the bisection method. The\n",
    "        initial interval [a_0,b_0] is given by [a,b]. If f(m_n) == 0 for some\n",
    "        midpoint m_n = (a_n + b_n)/2, then the function returns this solution.\n",
    "        If all signs of values f(a_n), f(b_n) and f(m_n) are the same at any\n",
    "        iteration, the bisection method fails and return None.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    >>> f = lambda x: x**2 - x - 1\n",
    "    >>> bisection(f,1,2,25)\n",
    "    1.618033990263939\n",
    "    >>> f = lambda x: (2*x - 1)*(x - 3)\n",
    "    >>> bisection(f,0,1,10)\n",
    "    0.5\n",
    "    '''\n",
    "    if f(a)*f(b) >= 0:\n",
    "        print(\"Bisection method fails.\")\n",
    "        return None\n",
    "    a_n = a\n",
    "    b_n = b\n",
    "    for n in range(1,N+1):\n",
    "        m_n = (a_n + b_n)/2\n",
    "        f_m_n = f(m_n)\n",
    "        if f(a_n)*f_m_n < 0:\n",
    "            a_n = a_n\n",
    "            b_n = m_n\n",
    "        elif f(b_n)*f_m_n < 0:\n",
    "            a_n = m_n\n",
    "            b_n = b_n\n",
    "        elif f_m_n == 0:\n",
    "            print(\"Found exact solution.\")\n",
    "            return m_n\n",
    "        else:\n",
    "            print(\"Bisection method fails.\")\n",
    "            return None\n",
    "    return (a_n + b_n)/2\n",
    "\n",
    "# Finds value of y for a given x using step size h \n",
    "# and initial value y0 at x0. \n",
    "def rungeKutta(x0, y0, x, h, arg, dydx):\n",
    "    \"\"\"\n",
    "    https://www.geeksforgeeks.org/runge-kutta-4th-order-method-solve-differential-equation/\n",
    "    \"\"\"\n",
    "    # Count number of iterations using step size or \n",
    "    # step height h \n",
    "    n = (int)((x - x0)/h)  \n",
    "    # Iterate for number of iterations \n",
    "    y = y0 \n",
    "    for i in range(1, n + 1): \n",
    "        \"Apply Runge Kutta Formulas to find next value of y\"\n",
    "        k1 = h * dydx(x0, y, arg) \n",
    "        k2 = h * dydx(x0 + 0.5 * h, y + 0.5 * k1, arg) \n",
    "        k3 = h * dydx(x0 + 0.5 * h, y + 0.5 * k2, arg) \n",
    "        k4 = h * dydx(x0 + h, y + k3, arg) \n",
    "  \n",
    "        # Update next value of y \n",
    "        y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4) \n",
    "  \n",
    "        # Update next value of x \n",
    "        x0 = x0 + h \n",
    "    return y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dadt(t, a,  a0half, OmegaR0=1e-5, OmegaM0=-1e-5+0.3, OmegaK0=0., OmegaLambda0=0.7):\n",
    "    \"\"\"\n",
    "     https://en.wikipedia.org/wiki/Friedmann_equations\n",
    "    \"\"\"\n",
    "    if OmegaK0 <0. and a> a0half:\n",
    "        return -sqrt(OmegaR0/(a*a)+OmegaM0/a+OmegaK0+OmegaLambda0*a*a)\n",
    "    else:\n",
    "        return sqrt(OmegaR0/(a*a)+OmegaM0/a+OmegaK0+OmegaLambda0*a*a)\n",
    "    \n",
    "def aMid(a, OmegaR0=1e-5,  OmegaM0=-1e-5+1.1, OmegaK0=-0.1, OmegaLambda0=0.):\n",
    "    \"\"\"\n",
    "    Exercise: \n",
    "        For OmegaK0<0 use this to find the maximum value a(t=tmax) can take. Use negative \n",
    "        square root in the integrand for t>tmax to.\n",
    "    \"\"\"\n",
    "    return OmegaR0/(a*a)+OmegaM0/a+OmegaK0+OmegaLambda0*a*a\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found exact solution.\n",
      "5.499954545492112\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "#x0 = 0\n",
    "#y = 1\n",
    "#x = 2\n",
    "#h = 0.2\n",
    "#print('The value of y at x is:', rungeKutta(x0, y, x, h))\n",
    "a0Mid = bisection(aMid, 1e-5,20.,1000)/2.\n",
    "print(a0Mid)\n",
    "print(aMid(a0Mid*2.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_init=1e-4#1e-7\n",
    "OmegaR0=1e-5\n",
    "a_init=sqrt(2.*sqrt(OmegaR0)*t_init)\n",
    "h=t_init/2.\n",
    "t_next=t_init*2.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0007952707287670507"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.002233795950365733"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rungeKutta(t_init, a_init, t_next, h, a0Mid, dadt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_vals = np.linspace(t_init,1.1,10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_vals = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "for t in t_vals:\n",
    "    a_vals.append(rungeKutta(t_init, a_init, t, h, a0Mid, dadt))\n",
    "    a_init = a_vals[-1]\n",
    "    t_init = t\n",
    "print(\"done\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1f880a7be0>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX+//HXSegQamhCQkIJEJpgAJUVRIqI2FcFdRUbrru4u7a1N0QE6+KKq4jYy9o30gVBbCABFEglJARCS+iQkDrn90fi9xcRzQgzc2du3s/HI4+dcjb3cx4zeXs599xzjLUWERFxlzCnCxAREd9TuIuIuJDCXUTEhRTuIiIupHAXEXEhhbuIiAsp3EVEXEjhLiLiQgp3EREXquXUgSMjI21MTIxThxcRCUmrV6/eba1tWV07x8I9JiaGpKQkpw4vIhKSjDE53rTTsIyIiAsp3EVEXEjhLiLiQgp3EREXUriLiLiQwl1ExIUU7iIiLqRwFxEJkMKSMh6fn8rWvYV+P5ZjNzGJiNQkS9PzuP+TDWzbf4T2zRrwp1M7+PV4CncRET/KP1TMpDkpfPbjdjq1bMj7N53GgNjmfj+uwl1ExA+stXyQlMtj81I5UlLOP4Z34eYzO1G3VnhAjl9tuBtjZgNjgDxrbc9jvG+A6cBooBAYb61d4+tCRURCRVb+Ye79ZD0rsvYyIKY5Uy7uSedWEQGtwZsz99eA54E3fuX9c4AulT8Dgf9U/q+ISI1SUubhpS838e+lmdStFcbjF/fi8oQowsJMwGupNtyttcuNMTG/0eQC4A1rrQVWGGOaGmPaWmt3+KhGEZGgtzpnL3d/tJ6NeYc5t3dbHhoTT6vG9Ryrxxdj7u2ArVWe51a+pnAXEdc7WFTKEwvSeGvFFk5qUo9XrklgWPfWTpcV2AuqxpgJwASA6OjoQB5aRMSnrLXMWbeDSXNS2HO4mOsGxXL7yDga1g2OeSq+qGIbEFXlefvK137BWjsTmAmQkJBgfXBsEZGAy9lTwAP/S2Z5Rj492zVm1tUJ9Ilq6nRZP+OLcE8EJhpj3qPiQuoBjbeLiBsVl5Uz88ssnl+aSe3wMB46L56rT4sh3IELptXxZirku8CZQKQxJhd4CKgNYK19EZhHxTTITCqmQl7rr2JFRJzy3aY93P/pejblF3Bur7Y8MCaeNk2cu2BaHW9my4yr5n0L/NVnFYmIBJHdh4uZMi+Vj9dsI6p5fV4d35+h3Vo5XVa1gmPkX0QkyHg8lv8mbWXq/DQKS8r469BOTBzahfp1AnOH6YlSuIuIHCVt50Hu+2QDq3P2MSC2OY9d2JMurQN7h+mJUriLiFQqLClj+uKNzPo6m8b1avHUpX24pF87KlZZCS0KdxERYHHKLh5KTGbb/iNcnhDF3ed0o1nDOk6XddwU7iJSo23bf4RHEpNZlLKLuNaN+ODPp9E/xv9L8vqbwl1EaqTisnJmfZXNv7/YCMBdo7pxwxmx1A53xwZ1CncRqXGWZ+TzUGIy2bsLGNWjDfeP6U77Zg2cLsunFO4iUmNs33+ER+ekMH/DTmJaNOC1a/tzZtfgn7N+PBTuIuJ6JWUeZn2dxb+XZGKx3DEyjhsHdwzYrkhOULiLiKt9vXE3DyZuICu/gJHxrXlgTDxRzd01BHMsCncRcaUdB44weU4qc9fvoEOLBiGzbICvKNxFxFVKyjzM/iab55ZspNxjuW1EHBMGd6RebfcOwRyLwl1EXOObzN08+L8NbMovYER8ax6sIUMwx6JwF5GQt/NAEZPnpjBn3Q6imzdg9vgEzurm/FZ3TlK4i0jIKinz8GrlEEyZx3Lr8DhuGlLzhmCOReEuIiFpWXoekz5LIWt3AcO6teKh83oQ3aJmDsEci8JdRELK5t0FTJ6bwuLUPGIjG9a4WTDeUriLSEgoKC5jxtJMZn2VTe1wwz3ndOPaQbHUqeWOtWB8TeEuIkHNWkvij9uZMi+VXQeLubhfO+4e1Y1WjYN3/9JgoHAXkaC1YdsBHvksmVWb99GrXRNeuPIUTunQzOmyQoLCXUSCzt6CEp5alM6732+heYM6TLukF5eeEkVYWOjtiOQUhbuIBI2ycg9vr9zC04vSKSgpZ/zpMfxjeBxN6td2urSQo3AXkaDw7abdPJKYQvquQwzq3IKHzutBXIhtSh1MFO4i4qht+48wZW7FAl/tm9Xnxav6cXaPNiG5KXUwUbiLiCOKSst56css/vNlJkCNXeDLXxTuIhJQ1lrmrt/B4/PS2Lb/COf2bsu9o7vTrml9p0tzFYW7iATMutz9TPoshaScfXRv25gnL+3N6Z0inS7LlRTuIuJ3uw4W8eTCdD5cnUtkozpMvbgXlyZEEa6pjX6jcBcRvykqLWfWV1m8sGwTZeWWm4Z0ZOLQzkTU09RGf1O4i4jPHT2ufnaP1tw7ujsdWjR0urQaw6twN8aMAqYD4cAsa+3Uo96PBl4Hmla2udtaO8/HtYpICFife4BJcyqWDNC4unOqDXdjTDgwAxgB5AKrjDGJ1tqUKs3uB9631v7HGBMPzANi/FCviASpn8bVP1qTS/MGdXj84l5cpnF1x3hz5j4AyLTWZgEYY94DLgCqhrsFGlc+bgJs92WRIhK8ikrLeeXrbGYszaS03MOEMzry17M601jj6o7yJtzbAVurPM8FBh7V5mFgkTHmFqAhMPxYv8gYMwGYABAdHf17axWRIHL0uPrI+Ipx9ZhIjasHA19dUB0HvGatfdoYcxrwpjGmp7XWU7WRtXYmMBMgISHB+ujYIhJgVcfVu7WJ4J0bBnJ6Z42rBxNvwn0bEFXlefvK16q6HhgFYK39zhhTD4gE8nxRpIgEh7yf5qtXjqtPuagXl/fXuHow8ibcVwFdjDGxVIT6WOCKo9psAYYBrxljugP1gHxfFioizjlSUs7LX2Xx4pebKC33cOMZHZmocfWgVm24W2vLjDETgYVUTHOcba1NNsZMApKstYnA7cDLxphbqbi4Ot5aq2EXkRDn8Vg+XruNpxams/NgEWf3aM0952hcPRR4NeZeOWd93lGvPVjlcQowyLeliYiTvt20m8fmppK8/SC92zdh+tiTGdixhdNliZd0h6qI/Mym/MM8Pi+Vxal5nNSkHv+6/GTO73OStrgLMQp3EQEq9i2dvjiDt1duoV7tcO48uyvX/yFW66uHKIW7SA1XVFrO699u5vmlmRQUlzFuQDS3jogjslFdp0uTE6BwF6mhrLXMWbeDaQvSyN13hKFdW3Lv6O500b6lrqBwF6mBVufsZfLcVNZu2U+3NhG8df1A/tBFNyG5icJdpAbZsqeQaQvSmLt+B60i6vLEJb255JT2ugnJhRTuIjXAgcJSnl+6kde/zSE8zPD3YV2YMLgjDesqAtxKn6yIi5WWe3hrRQ7Tl2zkwJFS/tivPbeP7EqbJvWcLk38TOEu4kLWWhal7GLq/DSydxcwqHML7h3dnR4nNXG6NAkQhbuIy6zPPcDkuSmszN5L51aNmD0+gaFdW2GMxtVrEoW7iEts3VvIkwvTSfxxOy0a1uHRC3syrn8UtcLDnC5NHKBwFwlx+wpKeH5pJm9+l0NYGEwc2pmbhnQkQis21mgKd5EQVVRazmvfbmZG5Z2ll54Sxa0j4nSxVACFu0jIKfdYPl27jacXpbP9QBFndWvFXaO60bWN7iyV/0/hLhJClmfk8/j8NFJ3VCzD+9RlfTi9k+4slV9SuIuEgOTtB5g6P42vNu4mqnl9nhvXlzG92moZXvlVCneRIJa7r5BnFmXwyQ/baFK/Ng+MieeqU6OpW0vL8MpvU7iLBKEDhaW8sCyTV7/dDMBNgztx85mdaFJfM2DEOwp3kSBSXFbOm9/l8O8vMjlYVMrFfdtz28g42jWt73RpEmIU7iJBwOOxfLZuO08uTCd33xGGxLXk7nO60b1tY6dLkxClcBdx2LeZu5kyP5UN2w7S46TGTL24t9ZWlxOmcBdxSNrOg0ydn8ay9HzaNa2vjajFpxTuIgG248ARnl6UwUdrcomoW4v7RnfnT6d10EbU4lMKd5EAOVhUyn+WbWL219lYCzee0ZG/nNmJpg3qOF2auJDCXcTPSso8vL0yh+eWbGRfYSkX9W3HbSPiiGrewOnSxMUU7iJ+Yq1l7vodPLEgnS17CxnUuQX3nNOdnu20YYb4n8JdxA9WZO3h8Xmp/Jh7gG5tInj9ugEM7hKpDTMkYBTuIj6UsesQ0+ansSQtj7ZN6vHUpX24qG87wjUDRgJM4S7iA7sOFvHs5xm8n7SVhnVqcdeoblw7KEYzYMQxXoW7MWYUMB0IB2ZZa6ceo81lwMOABX601l7hwzpFgtKholJmLs/i5a+yKPdYxp8ey8SzOtO8oWbAiLOqDXdjTDgwAxgB5AKrjDGJ1tqUKm26APcAg6y1+4wxrfxVsEgwKC338O73W5i+eCN7Cko4r89J3DmyK9EtNANGgoM3Z+4DgExrbRaAMeY94AIgpUqbG4EZ1tp9ANbaPF8XKhIMrLUs2LCTJxamk727gIGxzZk9ujt9opo6XZrIz3gT7u2ArVWe5wIDj2oTB2CM+YaKoZuHrbULfFKhSJBYtXkvU+alsnbLfrq0asTs8QkM7dpKM2AkKPnqgmotoAtwJtAeWG6M6WWt3V+1kTFmAjABIDo62keHFvGvzLzDPLEgjUUpu2jduC7TLunFJf3aUys8zOnSRH6VN+G+DYiq8rx95WtV5QIrrbWlQLYxJoOKsF9VtZG1diYwEyAhIcEeb9EigZB3qIjpizfy3qqt1K8dzh0j47juD7E0qKNJZhL8vPmWrgK6GGNiqQj1scDRM2E+BcYBrxpjIqkYpsnyZaEigVJQXMbLX2Uxc3kWJWUerhwYzd+GdSGyUV2nSxPxWrXhbq0tM8ZMBBZSMZ4+21qbbIyZBCRZaxMr3xtpjEkByoE7rbV7/Fm4iK+VlXv4b9JWnv18I7sPF3Nur7bccXZXYiMbOl2ayO9mrHVmdCQhIcEmJSU5cmyRqqy1LErZxbQFaWTlF9A/phn3jO5Ov+hmTpcm8gvGmNXW2oTq2mnwUGq01Tn7eHxeKkk5++jUsiEvX53A8O6aASOhT+EuNVL27gKeWJDG/A07aRlRlykX9eKyBM2AEfdQuEuNsvtwMc8t2cg7K7dQp1YYtw6P44YzYmlYV38K4i76RkuNUFhSxitfZfPil5soKvMwbkAUfx8WR8sIzYARd1K4i6uVlXv4YHUuz36eQd6hYs7u0Zp/jupGp5aNnC5NxK8U7uJK1lqWpOYxdUEamXmH6RfdlBeu7EdCTHOnSxMJCIW7uM4PW/fz+LxUVmbvpWNkQ168qh9n92ijGTBSoyjcxTVy9hTwxMJ05q7bQWSjOjx6YU/G9o+itmbASA2kcJeQt7eghOeWbOTtlTnUCgvjb8O6MGFwRxppBozUYPr2S8g6UlLO7G+yeXHZJgpKyri8fzS3Du9Cq8b1nC5NxHEKdwk55R7LR2tyeWZRBjsPFjG8e2vuPqcrnVtFOF2aSNBQuEvIsNayLCOfqfPSSN91iD5RTZk+9mQGdmzhdGkiQUfhLiEhZftBpsxL5evM3XRo0YAXruzHOT01A0bk1yjcJajtPFDE04vS+XBNLk3q1+bBMfFcdWoH6tTSDBiR36Jwl6BUUFzGS8uzeHl5FuUey41ndOSvZ3amSYPaTpcmEhIU7hJUyj2WD5K28vTnGeQfKmZM77b88+xuRLdo4HRpIiFF4S5B48uMfKbMTSV91yFO6dCMl/50ijbMEDlOCndxXNrOg0yZl8byjHyim+tiqYgvKNzFMXkHi3jm8wzeT9pKRL3a3H9ud/50Wgfq1gp3ujSRkKdwl4ArLClj5vIsZi7PorTcw7WDYrnlrM40bVDH6dJEXEPhLgFT7rF8tDqXpxalk3eomNG92nDXqG50aNHQ6dJEXEfhLgHx1cZ8HpubStrOQ/SNbsp/rurHKR20trqIvyjcxa8y8w7z2NwUlqbnE9W8Ps9f0Zdze7XVxVIRP1O4i1/sLyzhX4s38taKHOrXDufe0d245vQYXSwVCRCFu/hUabmHd1Zu4dnFGRw8Usq4AdHcNiKOFo20EbVIICncxWeWpecxeW4qmXmHGdS5BQ+Miadbm8ZOlyVSIync5YRl5h1i8txUlqXnExvZkFlXJzCseyuNq4s4SOEux+2ncfU3V+TQoE4495/bnatPi9GKjSJBQOEuv1tpuYe3V+Tw7OKNHCoq5YqB0dw6XOPqIsHEq3A3xowCpgPhwCxr7dRfaXcJ8CHQ31qb5LMqJWgsTc9j8pwUNuUX8IfOkTwwJp6ubbS9nUiwqTbcjTHhwAxgBJALrDLGJFprU45qFwH8HVjpj0LFWZl5h3h0TipfZlSMq79yTQJnddO4ukiw8ubMfQCQaa3NAjDGvAdcAKQc1e5RYBpwp08rFEcdLCrlucUbee3bzdTXuLpIyPAm3NsBW6s8zwUGVm1gjOkHRFlr5xpjFO4u4PFYPlqTy7QF6ewpKGZs/2juGKlxdZFQccIXVI0xYcAzwHgv2k4AJgBER0ef6KHFT37cup+HEpP5Yet++kU35dXx/enVvonTZYnI7+BNuG8Doqo8b1/52k8igJ7Assrx1zZAojHm/KMvqlprZwIzARISEuwJ1C1+sPtwMU8sSOP9pFxaRtTlmcv6cOHJ7QgL07i6SKjxJtxXAV2MMbFUhPpY4Iqf3rTWHgAif3pujFkG3KHZMqGjtNzDm9/l8OziDI6UlDNhcEduOaszEfW0GbVIqKo23K21ZcaYicBCKqZCzrbWJhtjJgFJ1tpEfxcp/vNt5m4e/iyZjF2HGRzXkgfHxNO5VSOnyxKRE+TVmLu1dh4w76jXHvyVtmeeeFnib7n7CnlsbirzN+wkunkDXr46geFaMkDENXSHag1TXFbOzC+zeH5pJsbAHSPjuOGMjtSrraV4RdxE4V6DfLUxnwf/l0z27gJG92rD/efGc1LT+k6XJSJ+oHCvAXYdLOLROSnMWbeDmBYNeP26AQyJa+l0WSLiRwp3Fysr9/D6dzk8+3kGJeUebh0ex01DNAQjUhMo3F1qdc5e7vtkA2k7D3Fm15Y8cn4POrRo6HRZIhIgCneX2VtQwrT5afw3aSttm9TjxatO4ewerTULRqSGUbi7hMdjeT9pK1MXpHG4qIybhnTkb2d1oWFdfcQiNZH+8l0gY9ch7vl4Patz9jEgtjmTL+xJXGutsS5SkyncQ1hRaTkzlmby4pebaFS3Fk9d2odL+rXTEIyIKNxD1YqsPdz78Xqydhdwcd923Hdudy3HKyL/R+EeYg4UlvL4/FTeW7WVqOb1eeO6AQzWnHUROYrCPURYa5mzbgePfJbCvsISbhrSkX8Mi6N+Hc1ZF5FfUriHgNx9hTz4v2S+SMujd/smvH5df3qcpM0zROTXKdyDmMdjeeO7zTyxMB2AB8bEM/70GMK1eYaIVEPhHqSydxdw14fr+H7zXobEteSxi3rSvlkDp8sSkRChcA8y5R7Lq99k89SidGqHh/HkH3vzx1Paa3qjiPwuCvcgsin/MHd+8CNrtuxnWLdWTLm4F60b13O6LBEJQQr3IFDuscz6KounP8+gfu1wnr28YmNqna2LyPFSuDts465D3PHhOn7cup+R8a2ZfFFPWkXobF1ETozC3SEej+WVr7N5clE6DeuE89y4vpzXu63O1kXEJxTuDsjdV8jt7//Iyuy9jIhvzeMX9yJSSweIiA8p3APIWstHa7bxSGIyFnjij725VDNhRMQPFO4BsreghHs/Xs+C5J0MiGnO05f1Iaq55q2LiH8o3APgi7Rd/PPD9Rw8Uso953TjhjM66i5TEfErhbsfHSkp59G5Kbyzcgvd2kTw5vUD6N62sdNliUgNoHD3k7SdB7nlnbVszDvMhMEduX1kHHVraQVHEQkMhbuPWWt5a+UWJs9JIaJeba23LiKOULj70P7CEu76aB0Lk3cxOK4lT1/ah5YRmuIoIoGncPeR77P38o/31pJ/uJj7Rnfn+j/EEqaLpiLikDBvGhljRhlj0o0xmcaYu4/x/m3GmBRjzDpjzBJjTAfflxqcPB7Lv5dsZOzM76hTK4yPbj6dGwd3VLCLiKOqPXM3xoQDM4ARQC6wyhiTaK1NqdJsLZBgrS00xtwMPAFc7o+Cg8m+ghJuff8HlqXnc8HJJ/HYRb1oVFf/GBIR53mTRAOATGttFoAx5j3gAuD/wt1au7RK+xXAVb4sMhity93PzW+tIf9QMZMv7MmVA6N1p6mIBA1vwr0dsLXK81xg4G+0vx6YfyJFBTNrLW+v3MKkz1JoGVGXD/58Gn2imjpdlojIz/h0DMEYcxWQAAz5lfcnABMAoqOjfXnogCgsKeO+TzbwydptDIlryb8uP5lmDes4XZaIyC94E+7bgKgqz9tXvvYzxpjhwH3AEGtt8bF+kbV2JjATICEhwf7uah20dW8hN76RRPquQ9w2Io6JQzvroqmIBC1vwn0V0MUYE0tFqI8FrqjawBjTF3gJGGWtzfN5lQ77btMe/vL2aso9lteuHcAQ3ZQkIkGu2nC31pYZYyYCC4FwYLa1NtkYMwlIstYmAk8CjYAPKi8qbrHWnu/HugPmzRU5PJKYTIcWDZh1TX9iIxs6XZKISLW8GnO31s4D5h312oNVHg/3cV2OKynz8PBnybyzcgtDu7Zk+ri+NK5X2+myRES8oknZx7CvoISb3lrN99l7+fOQTtx5dlct0SsiIUXhfpScPQWMf3UV2/Yf4V+Xn8yFfds5XZKIyO+mcK9i7ZZ93PB6EuXW8s4NA0mIae50SSIix0XhXmlh8k7+/t5aWkXU47Vr+9OxZSOnSxIROW4Kd+C1b7J5ZE4Kfdo3ZdY1CUQ20jK9IhLaanS4W2t59vMMnvsik5HxrZk+ti/162i3JBEJfTU23D0ey6Q5Kbz27WYuT4hiysW9NCNGRFyjRoZ7WbmHf360jo/XbOOGP8Ry37ndtaKjiLhKjQv34rJy/vbuWhYm7+L2EXFMPKuzgl1EXKdGhXtJmYe/vr2Wxam7eOi8eK4dFOt0SSIiflFjwr203MMt765hceouHr2gB386LcbpkkRE/MarPVRDXVm5h7+/VzEU89B58Qp2EXE914e7x2O5/YMfmbd+J/ef211DMSJSI7g63K21PDo3hf/9sJ07z+7KDWd0dLokEZGAcHW4v/hlFq9+s5nrBsXylzM7OV2OiEjAuDbcP1ydy7QFaZzX5yTu1zx2EalhXBnuK7P2cPdH6xjUuQVPXdpbe52KSI3junDP3VfIzW+vIbp5A1648hTq1tJaMSJS87gq3AtLyrjxjdWUlnt4+ZoEmtTXtngiUjO5Ktzv+2QD6TsP8ty4vnTSeuwiUoO5Jtw/XpPLJ2u38bdhXRjatZXT5YiIOMoV4Z69u4AHPt3AgNjm3HJWF6fLERFxXMiHe7nHctv7P1C7VhjTx56sNdlFRHDBwmFvfreZtVv286/LT6Ztk/pOlyMiEhRC+sx9+/4jPLkwncFxLbng5JOcLkdEJGiEdLg/uTCdMo/lsQt76g5UEZEqQjbcN2w7wKc/bOPaQbFENW/gdDkiIkElZMP92c8zaFyvNjdrQTARkV8IyXDPyj/MkrQ8xp8eo7tQRUSOwatwN8aMMsakG2MyjTF3H+P9usaY/1a+v9IYE+PrQqt647sc6oSHcdWpHfx5GBGRkFVtuBtjwoEZwDlAPDDOGBN/VLPrgX3W2s7As8A0Xxf6E4/HMm/9Ds7q1oqWEXX9dRgRkZDmzZn7ACDTWptlrS0B3gMuOKrNBcDrlY8/BIYZP01fSd5+kLxDxYzq2cYfv15ExBW8Cfd2wNYqz3MrXztmG2ttGXAAaOGLAo+2YfsBAPpFN/PHrxcRcYWAXlA1xkwwxiQZY5Ly8/OP63e0aFiHEfGtad9Md6OKiPwab5Yf2AZEVXnevvK1Y7XJNcbUApoAe47+RdbamcBMgISEBHs8BY/s0YaRPTQkIyLyW7w5c18FdDHGxBpj6gBjgcSj2iQC11Q+/iPwhbX2uMJbREROXLVn7tbaMmPMRGAhEA7MttYmG2MmAUnW2kTgFeBNY0wmsJeK/wCIiIhDvFoV0lo7D5h31GsPVnlcBFzq29JEROR4heQdqiIi8tsU7iIiLqRwFxFxIYW7iIgLKdxFRFzIODUd3RiTD+Qc5/89Etjtw3KCjZv7p76FLjf3L5T61sFa27K6Ro6F+4kwxiRZaxOcrsNf3Nw/9S10ubl/buybhmVERFxI4S4i4kKhGu4znS7Az9zcP/UtdLm5f67rW0iOuYuIyG8L1TN3ERH5DUEd7sG2MbcvedG324wxKcaYdcaYJcaYkNoNvLr+VWl3iTHGGmNCZqaCN30zxlxW+fklG2PeCXSNx8uL72W0MWapMWZt5XdztBN1Hg9jzGxjTJ4xZsOvvG+MMc9V9n2dMaZfoGv0KWttUP5QsbzwJqAjUAf4EYg/qs1fgBcrH48F/ut03T7s21CgQeXjm0Olb972r7JdBLAcWAEkOF23Dz+7LsBaoFnl81ZO1+3Dvs0Ebq58HA9sdrru39G/wUA/YMOvvD8amA8Y4FRgpdM1n8hPMJ+5B9XG3D5Wbd+stUuttYWVT1dQsQNWqPDmswN4FJgGFAWyuBPkTd9uBGZYa/cBWGvzAlzj8fKmbxZoXPm4CbA9gPWdEGvtcir2m/g1FwBv2AorgKbGmLaBqc73gjncg2pjbh/zpm9VXU/FGUWoqLZ/lf/kjbLWzg1kYT7gzWcXB8QZY74xxqwwxowKWHUnxpu+PQxcZYzJpWKPh1sCU1pA/N6/y6Dm1WYd4hxjzFVAAjDE6Vp8xRgTBjwDjHe4FH+pRcXQzJlU/ItruTGml7V2v6NV+cY44DVr7dPGmNOo2IGtp7XW43Rh8nPBfOb+ezbm5rc25g5C3vQNY8xw4D7gfGttcYBq84Xq+hcB9ASWGWM2UzG+mRgiF1W9+exygURrbam1NhvIoCLsg503fbseeB/AWvsdUI+KdVncwKu/y1ARzOHu5o25q+2bMaYv8BIVwR4qY7Y/+c3+WWsPWGsjrbV5wSjsAAAA8klEQVQx1toYKq4pnG+tTXKm3N/Fm+/lp1SctWOMiaRimCYrkEUeJ2/6tgUYBmCM6U5FuOcHtEr/SQSurpw1cypwwFq7w+mijpvTV3R/64eKq9cZVFzBv6/ytUlUBAFUfLE+ADKB74GOTtfsw74tBnYBP1T+JDpdsy/7d1TbZYTIbBkvPztDxbBTCrAeGOt0zT7sWzzwDRUzaX4ARjpd8+/o27vADqCUin9dXQ/8Gfhzlc9tRmXf14fSd/JYP7pDVUTEhYJ5WEZERI6Twl1ExIUU7iIiLqRwFxFxIYW7iIgLKdxFRFxI4S4i4kIKdxERF/p/iTU2V2u4WXUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(t_vals, a_vals)\n",
    "#ax.set_yscale('log')\n",
    "#ax.set_xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}