{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Week 10 Day 3: ODE problems\n", "\n", "## Objectives\n", "\n", "* Discuss debuggers\n", "* Solve a BVP\n", "* Look at some other ODE problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Debugging\n", "\n", "Using a debugger has several parts:\n", "\n", "* Starting the debugger\n", " * Classic way: `import pdb; pdb.set_trace()` (This is the standard library debugger)\n", " * Python 3.7 way: `breakpoint()` (This will start the current debugger)\n", " * IPython offers `%debug` to jump in *after* an error is raised! (like `pdb.pm()`)\n", " * You can also start at the beginning (see [docs](https://docs.python.org/3.7/library/pdb.html))\n", "* Controlling the debugger\n", " * `pyb` (and `ipydb`) are sort of based on classic debugging tools - the syntax is a bit odd, but will be useful in other languages.\n", " * A common term is the \"stack\", short in this case for the \"call-stack\", the series of functions that call each other to get to the current point.\n", " * Note that the \"stack\" is also a term for memory locations in compiled programs - and is unrelated to this usage." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| short | name | action | Names in other debuggers |\n", "|-------|---------------|--------|--------------------------|\n", "| `h` | `help` | Print out help | |\n", "| `p` | `print` | Print out the value of a variable or expression | |\n", "| `u` | `up` | Go up one in the stack| |\n", "| `d` | `down` | Go down one in the stack | |\n", "| `w` | `where` | Show you where you are | Stack trace |\n", "| `s` | `step` | Move forward one computation | Step in |\n", "| `n` | `next` | More forward one line | Step over |\n", "| `c` | `continue` | Continue on till next stop | |\n", "| `b` | `breakpoint` | Set a breakpoint somehere | |\n", "| `q` | `quit` | Quit | |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## A program to debug\n", "def simple_function(div):\n", " simple_x = 2\n", " simple_y = div\n", "\n", " # import pdb; pdb.set_trace()\n", "\n", " return fancy_function(simple_x, simple_y)\n", "\n", "\n", "def fancy_function(x, y):\n", " r = x / y\n", " return r" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# simple_function(0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %debug" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# simple_function(2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import pdb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# pdb.run('simple_function(0)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Other tools: Tracing\n", "\n", "Some IDEs, like PyCharm, may offer enhanced and partially graphical debuggers. And, in Python 3.7, the new built-in `breakpoint()` can call the fancy debugger instead of the builtin one!\n", "\n", "Let's look at something similar: Tracing. If you have a complex piece of Python code and you want to have an idea of what the control flow looks like, you can trace it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%writefile simple.py\n", "\n", "## A program to debug\n", "def simple_function(div):\n", " simple_x = 2\n", " simple_y = div\n", "\n", " return fancy_function(simple_x, simple_y)\n", "\n", "\n", "def fancy_function(x, y):\n", " r = x / y\n", " return r\n", "\n", "\n", "def main():\n", " print(simple_function(2))\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!python simple.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!python -m trace --trace simple.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!python -m trace --listfuncs simple.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary value problems (BPV)\n", "\n", "Let's look at a square well:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.integrate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align}\n", "\\psi'' + \\left( \\frac{2 m}{\\hbar^2} V_0 - \\kappa^2 \\right) \\psi = 0,& \\quad \\textrm{for }|x|\\le a \\\\\n", "\\psi'' - \\kappa^2 \\psi = 0,&\n", "\\quad \\textrm{for }|x| > a\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, we write this as a set of first order equations:\n", "\n", "$$\n", "u = \\left(\\begin{matrix}\\psi' \\\\ \\psi\\end{matrix}\\right)\n", "$$\n", "\n", "$$\n", "u' = \\left(\\begin{matrix}\\psi'' \\\\ \\psi'\\end{matrix}\\right)\n", "= \\left(\\begin{matrix}\\left( \\frac{2 m}{\\hbar^2} V_0 - \\kappa^2 \\right) \\psi \\\\ \\psi'\\end{matrix}\\right)\n", "$$\n", "\n", "where we set $V_0$ to $0$ when $|x|$ is larger than $a$. Boundary conditions:\n", "\n", "$$\n", "\\psi(x_\\mathrm{max}) = e^{-\\kappa x_\\mathrm{max}} \\\\\n", "\\psi(-x_\\mathrm{max}) = e^{-\\kappa x_\\mathrm{max}}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Double click to see broken code for `solve_bvp` (in a live notebook, not a viewer).\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased)\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2017. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# QuantumEigen.py: Finds E and psi via rk4 + bisection\n", "\n", "# mass/((hbar*c)**2)= 940MeV/(197.33MeV-fm)**2 =0.4829\n", "# well width=20.0 fm, depth 10 MeV, Wave function not normalized.\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "eps = 1e-3 # Precision\n", "n_steps = 501\n", "E = -17 # E guess\n", "h = 0.04\n", "count_max = 100\n", "Emax = 1.1 * E # E limits\n", "Emin = E / 1.1\n", "\n", "\n", "def f(x, y, E):\n", " return np.array([y[1], -0.4829 * (E - V(x)) * y[0]])\n", "\n", "\n", "def V(x):\n", " # Well depth\n", " return -16 if abs(x) < 10 else 0\n", "\n", "\n", "def rk4(t, y, h, Neqs, E):\n", " ydumb = np.zeros((Neqs), float)\n", " k1 = np.zeros((Neqs), float)\n", " k2 = np.zeros((Neqs), float)\n", " k3 = np.zeros((Neqs), float)\n", " k4 = np.zeros((Neqs), float)\n", "\n", " F = f(t, y, E)\n", " k1 = h * F\n", " ydumb = y + k1 / 2\n", "\n", " F = f(t + h / 2, ydumb, E)\n", " k2 = h * F\n", " ydumb = y + k2 / 2\n", "\n", " F = f(t + h / 2, ydumb, E)\n", " k3 = h * F\n", " ydumb = y + k3\n", "\n", " F = f(t + h, ydumb, E)\n", " k4 = h * F\n", "\n", " y += (k1 + 2 * (k2 + k3) + k4) / 6\n", "\n", "\n", "def diff(E, h):\n", " i_match = n_steps // 3 # Matching radius\n", " nL = i_match + 1\n", "\n", " y = np.array([1e-15, 1e-15 * np.sqrt(-E * 0.4829)]) # Initial left wf\n", "\n", " for ix in range(0, nL + 1):\n", " x = h * (ix - n_steps / 2)\n", " rk4(x, y, h, 2, E)\n", "\n", " left = y[1] / y[0] # Log derivative\n", "\n", " y[0] = 1e-15\n", " # For even; reverse for odd\n", " y[1] = -y[0] * np.sqrt(-E * 0.4829) # Initialize R wf\n", "\n", " for ix in range(n_steps, nL + 1, -1):\n", " x = h * (ix + 1 - n_steps / 2)\n", " rk4(x, y, -h, 2, E)\n", "\n", " right = y[1] / y[0] # Log derivative\n", " return (left - right) / (left + right)\n", "\n", "\n", "def plot(E, h):\n", " # Repeat integrations for plot\n", " x = 0.0\n", " n_steps = 1501 # integration steps\n", "\n", " y = np.zeros((2), float)\n", " yL = np.zeros((2, 505), float)\n", " i_match = 500 # Matching point\n", " nL = i_match + 1\n", "\n", " y[0] = 1e-40 # Initial left wf\n", " y[1] = -np.sqrt(-E * 0.4829) * y[0]\n", "\n", " for ix in range(nL + 1):\n", " yL[0][ix] = y[0]\n", " yL[1][ix] = y[1]\n", " x = h * (ix - n_steps / 2)\n", " rk4(x, y, h, 2, E)\n", "\n", " y[0] = -1.0e-15 # For even; reverse for odd\n", " y[1] = -np.sqrt(-E * 0.4829) * y[0]\n", " j = 0\n", "\n", " Rwf_x = np.zeros(n_steps - 3 - nL)\n", " Rwf_y = np.zeros(n_steps - 3 - nL)\n", " for ix in range(n_steps - 1, nL + 2, -1): # right WF\n", " x = h * (ix + 1 - n_steps / 2) # Integrate in\n", " rk4(x, y, -h, 2, E)\n", " Rwf_x[j] = 2 * (ix + 1 - n_steps / 2)\n", " Rwf_y[j] = y[0] * 35e-9\n", " j += 1\n", " x = x - h\n", " normL = y[0] / yL[0][nL]\n", " j = 0\n", "\n", " # Renormalize L wf & derivative\n", " Lwf_x = np.zeros(nL + 1)\n", " Lwf_y = np.zeros(nL + 1)\n", " for ix in range(nL + 1):\n", " x = h * (ix - n_steps / 2 + 1)\n", " y[0] = yL[0][ix] * normL\n", " y[1] = yL[1][ix] * normL\n", " Lwf_x[j] = 2 * (ix - n_steps / 2 + 1)\n", " Lwf_y[j] = y[0] * 35e-9 # Factor for scale\n", " j += 1\n", "\n", " print(E)\n", " plt.plot(Rwf_x, Rwf_y)\n", " plt.plot(Lwf_x, Lwf_y)\n", " plt.show()\n", "\n", "\n", "for count in range(0, count_max + 1):\n", " # Iteration loop\n", " E = (Emax + Emin) / 2 # Divide E range\n", " Diff = diff(E, h)\n", "\n", " if diff(Emax, h) * Diff > 0:\n", " Emax = E # Bisection algor\n", " else:\n", " Emin = E\n", "\n", " if abs(Diff) < eps:\n", " break\n", "\n", " plot(E, h)\n", "\n", "\n", "print(\"Final eigenvalue E = \", E)\n", "print(\"iterations, max = \", count)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 2: Projectile motion with drag\n", "\n", "We can select a power $n=1$ or $2$. The friction coefficient is $k$.\n", "\n", "$$\n", "\\ddot{x} = - k \\dot{x}^n \\frac{\\dot{x}}{v} \\\\\n", "\\ddot{y} = - g - k \\dot{y}^n \\frac{\\dot{y}}{v} \\\\\n", "v = \\sqrt{\\dot{x}^2 + \\dot{y}^2}\n", "$$\n", "\n", "We can rewrite it:\n", "\n", "$$\n", "u = \\left(\\begin{matrix}\n", " x \\\\\n", " \\dot{x} \\\\\n", " y \\\\\n", " \\dot{y}\n", "\\end{matrix}\\right)\n", "$$\n", "\n", "Then:\n", "\n", "$$\n", "\\dot{u} = \\left(\\begin{matrix}\n", " \\dot{x} \\\\\n", " - k \\dot{x}^n \\frac{\\dot{x}}{v} \\\\\n", " \\dot{y} \\\\\n", " - g - k \\dot{y}^n \\frac{\\dot{y}}{v}\n", "\\end{matrix}\\right)\n", "=\n", "\\left(\\begin{matrix}\n", " u_1 \\\\\n", " - k u_1^n \\frac{u_1}{\\sqrt{u_1^2 + u_3^2}} \\\\\n", " u_3 \\\\\n", " - g - k u_3^n \\frac{u_3}{\\sqrt{u_1^2 + u_3^2}}\n", "\\end{matrix}\\right)\n", "$$\n", "\n", "We can set the initial positions to 0 and give it a starting velocity, and compare no friction to friction versions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.integrate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "v0 = 22\n", "angle = 34\n", "g = 9.8\n", "k = 0.8\n", "n = 1\n", "\n", "v0x = v0 * np.cos(np.radians(angle))\n", "v0y = v0 * np.sin(np.radians(angle))\n", "\n", "t_eval = np.linspace(0, 2)\n", "analytic_x = v0x * t_eval\n", "analytic_y = v0y * t_eval - 0.5 * g * t_eval ** 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def f(t, u):\n", " v = np.sqrt(u[1] ** 2 + u[3] ** 2)\n", " return np.stack(\n", " [u[1], -k * u[1] ** n * u[1] / v, u[3], -g - k * u[3] ** n * u[3] / v]\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = scipy.integrate.solve_ivp(f, [0, 2], [0, v0x, 0, v0y], t_eval=t_eval)\n", "print(res.message)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8, 3))\n", "plt.plot(res.y[0], res.y[2])\n", "plt.plot(analytic_x, analytic_y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased)\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2017. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# ProjectileAir.py: Order dt^2 projectile trajectory + drag\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "v0 = 22\n", "angle = 34.0\n", "g = 9.8\n", "kf = 0.8\n", "N = 20\n", "\n", "v0x = v0 * np.cos(angle * np.pi / 180)\n", "v0y = v0 * np.sin(angle * np.pi / 180)\n", "T = 2 * v0y / g\n", "H = v0y ** 2 / 2 / g\n", "R = 2 * v0x * v0y / g\n", "\n", "print(\"No Drag T =\", T, \", H =\", H, \", R =\", R)\n", "\n", "\n", "def plotNumeric(k):\n", " vx = v0 * np.cos(angle * np.pi / 180)\n", " vy = v0 * np.sin(angle * np.pi / 180)\n", " x = 0.0\n", " y = 0.0\n", " dt = vy / g / N * 1.5\n", "\n", " print(\"\\n With Friction \")\n", " print(\" x y\")\n", "\n", " xy = np.empty((2, N))\n", " xy[:, 0] = 0, 0\n", "\n", " for i in range(N - 1):\n", " vx = vx - k * vx * dt\n", " vy = vy - g * dt - k * vy * dt\n", " x = x + vx * dt\n", " y = y + vy * dt\n", " xy[:, i + 1] = x, y\n", " print(\" %13.10f %13.10f \" % (x, y))\n", "\n", " plt.plot(*xy)\n", "\n", "\n", "def plotAnalytic():\n", " v0x = v0 * np.cos(angle * np.pi / 180)\n", " v0y = v0 * np.sin(angle * np.pi / 180)\n", " dt = v0y / g / N * 1.8\n", " print(\"\\n No Friction \")\n", " print(\" x y\")\n", "\n", " xy = np.empty((2, N))\n", " for i in range(N):\n", " t = i * dt\n", " x = v0x * t\n", " y = v0y * t - g * t * t / 2\n", " xy[:, i] = x, y\n", " print(\" %13.10f %13.10f\" % (x, y))\n", "\n", " plt.plot(*xy)\n", "\n", "\n", "plotNumeric(kf)\n", "plotAnalytic()" ] } ], "metadata": { "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.8.2" } }, "nbformat": 4, "nbformat_minor": 2 }