variational-hw7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational principle\n",
    "Code implementing the variational principle for a 1D potential.\n",
    "To try different things, change:\n",
    "- the potential function V\n",
    "- the wavefunction psi\n",
    "- the initial guess of the wavefunction parameters aguess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This file has the solution to the last part of HW7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:29:18.220188Z",
     "start_time": "2018-10-19T18:29:18.188940Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numdifftools as nd\n",
    "from scipy.optimize import minimize\n",
    "from scipy.integrate import quad\n",
    "#from scipy.misc import derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:29:18.251437Z",
     "start_time": "2018-10-19T18:29:18.220188Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# some constants.  Changing mass will also change the energy scale.\n",
    "#hbar = 1.055e-34 # Planck's constant (J s)\n",
    "#mass = 1e-30 #  electron mass (kg)\n",
    "#L = 1e-9 # well half-width (m)\n",
    "#V0 = 1.6e-19 # well depth (J)\n",
    "# rescaled units:\n",
    "hbar =1.0\n",
    "mass = 1.0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the problem we're solving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:43:14.569183Z",
     "start_time": "2018-10-19T18:43:14.553560Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define the potential energy function\n",
    "\n",
    "# some constants for the potential\n",
    "L = 1.0\n",
    "V0 = 18*hbar**2/(2*mass*L**2)\n",
    "\n",
    "# domain of the potential\n",
    "xmin = -10*L; xmax = 10*L #limits of integration, may be np.inf (infinity)\n",
    "\n",
    "# The potential energy function itself\n",
    "def V(x):\n",
    "    # this example is the finite square well\n",
    "    # make sure this function is \"vectorizable\" and does not have any division by zero\n",
    "    return np.piecewise(x, [np.abs(x)<= L, np.abs(x)>L],[0.0,V0])\n",
    "\n",
    "bad_points = [-L,L] # list of x values where discontinuities happen in V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:43:45.224469Z",
     "start_time": "2018-10-19T18:43:45.208844Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Define the trial ground state wavefunction\n",
    "def psi(x,a):\n",
    "    '''Wavefunction as a function of position x.\n",
    "    The array \"a\" stores any parameters.'''\n",
    "    # Note, normalization is not necessary\n",
    "    #return np.exp(-x**2*a[0]) # gaussian\n",
    "    return np.piecewise(x,[np.abs(x*a[0])<1.0,np.abs(x*a[0])>=1.0],[1.0+np.cos(a[0]*np.pi*x),0.0])# cosine piece"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:00.692735Z",
     "start_time": "2018-10-19T18:44:00.677111Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "a_guess = [0.95/L] # initial guess of the arguments that minimize the wavefunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the calculations..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:01.536461Z",
     "start_time": "2018-10-19T18:44:01.520834Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def Hpsi(psi,x,a):\n",
    "    '''The Hamiltonian acting on psi at point x'''\n",
    "    return -0.5*hbar**2/mass *nd.Derivative(lambda x: psi(x,a), n=2)(x) + V(x)*psi(x,a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:01.802076Z",
     "start_time": "2018-10-19T18:44:01.770828Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def EVH(a, psi):\n",
    "    '''Expectation value of the Hamiltonian for wavefunction psi that has the arguments a\n",
    "    psi is a function of x and a\n",
    "    a is a list of parameters for the wavefunction'''\n",
    "    # Calculate normalization of psi\n",
    "    norm = quad(lambda x :np.abs(psi(x,a))**2, xmin, xmax, points=bad_points)[0]\n",
    "    ev = quad(lambda x :np.conj(psi(x,a))*Hpsi(psi, x,a), xmin, xmax,points=bad_points)[0]\n",
    "    return ev/norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:02.317686Z",
     "start_time": "2018-10-19T18:44:02.223936Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy of guess = 1.4845628466660217\n"
     ]
    }
   ],
   "source": [
    "# Calculate the expectation value of energy for our initial guess\n",
    "initEV = EVH(a_guess,psi)\n",
    "print(\"Initial energy of guess = {}\".format(initEV))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:02.317686Z",
     "start_time": "2018-10-19T18:44:02.223936Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[10], line 4\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# minimize the expectation value of energy by changing a_guess\u001b[39;00m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;66;03m# `output` has two pieces\u001b[39;00m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;66;03m# `output.x` is the list of parameters (what I called a above)\u001b[39;00m\n\u001b[0;32m----> 4\u001b[0m output \u001b[38;5;241m=\u001b[39m \u001b[43mminimize\u001b[49m\u001b[43m(\u001b[49m\u001b[43mEVH\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43ma_guess\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_minimize.py:691\u001b[0m, in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    689\u001b[0m     res \u001b[38;5;241m=\u001b[39m _minimize_cg(fun, x0, args, jac, callback, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n\u001b[1;32m    690\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbfgs\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m--> 691\u001b[0m     res \u001b[38;5;241m=\u001b[39m \u001b[43m_minimize_bfgs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfun\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mjac\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcallback\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    692\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnewton-cg\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m    693\u001b[0m     res \u001b[38;5;241m=\u001b[39m _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,\n\u001b[1;32m    694\u001b[0m                              \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_optimize.py:1388\u001b[0m, in \u001b[0;36m_minimize_bfgs\u001b[0;34m(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, xrtol, **unknown_options)\u001b[0m\n\u001b[1;32m   1385\u001b[0m pk \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m-\u001b[39mnp\u001b[38;5;241m.\u001b[39mdot(Hk, gfk)\n\u001b[1;32m   1386\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m   1387\u001b[0m     alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 \u001b[38;5;241m=\u001b[39m \\\n\u001b[0;32m-> 1388\u001b[0m              \u001b[43m_line_search_wolfe12\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmyfprime\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgfk\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1389\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43mold_fval\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mold_old_fval\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mamin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-100\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mamax\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e100\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1390\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m _LineSearchError:\n\u001b[1;32m   1391\u001b[0m     \u001b[38;5;66;03m# Line search failed to find a better solution.\u001b[39;00m\n\u001b[1;32m   1392\u001b[0m     warnflag \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_optimize.py:1160\u001b[0m, in \u001b[0;36m_line_search_wolfe12\u001b[0;34m(f, fprime, xk, pk, gfk, old_fval, old_old_fval, **kwargs)\u001b[0m\n\u001b[1;32m   1146\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1147\u001b[0m \u001b[38;5;124;03mSame as line_search_wolfe1, but fall back to line_search_wolfe2 if\u001b[39;00m\n\u001b[1;32m   1148\u001b[0m \u001b[38;5;124;03msuitable step length is not found, and raise an exception if a\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1155\u001b[0m \n\u001b[1;32m   1156\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1158\u001b[0m extra_condition \u001b[38;5;241m=\u001b[39m kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mextra_condition\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[0;32m-> 1160\u001b[0m ret \u001b[38;5;241m=\u001b[39m \u001b[43mline_search_wolfe1\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfprime\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgfk\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1161\u001b[0m \u001b[43m                         \u001b[49m\u001b[43mold_fval\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mold_old_fval\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1162\u001b[0m \u001b[43m                         \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1164\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ret[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m extra_condition \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   1165\u001b[0m     xp1 \u001b[38;5;241m=\u001b[39m xk \u001b[38;5;241m+\u001b[39m ret[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m*\u001b[39m pk\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_linesearch.py:84\u001b[0m, in \u001b[0;36mline_search_wolfe1\u001b[0;34m(f, fprime, xk, pk, gfk, old_fval, old_old_fval, args, c1, c2, amax, amin, xtol)\u001b[0m\n\u001b[1;32m     80\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39mdot(gval[\u001b[38;5;241m0\u001b[39m], pk)\n\u001b[1;32m     82\u001b[0m derphi0 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mdot(gfk, pk)\n\u001b[0;32m---> 84\u001b[0m stp, fval, old_fval \u001b[38;5;241m=\u001b[39m \u001b[43mscalar_search_wolfe1\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m     85\u001b[0m \u001b[43m        \u001b[49m\u001b[43mphi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mderphi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mold_fval\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mold_old_fval\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mderphi0\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m     86\u001b[0m \u001b[43m        \u001b[49m\u001b[43mc1\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mc1\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mc2\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mc2\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mamax\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mamax\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mamin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mamin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mxtol\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     88\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m stp, fc[\u001b[38;5;241m0\u001b[39m], gc[\u001b[38;5;241m0\u001b[39m], fval, old_fval, gval[\u001b[38;5;241m0\u001b[39m]\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_linesearch.py:161\u001b[0m, in \u001b[0;36mscalar_search_wolfe1\u001b[0;34m(phi, derphi, phi0, old_phi0, derphi0, c1, c2, amax, amin, xtol)\u001b[0m\n\u001b[1;32m    159\u001b[0m     alpha1 \u001b[38;5;241m=\u001b[39m stp\n\u001b[1;32m    160\u001b[0m     phi1 \u001b[38;5;241m=\u001b[39m phi(stp)\n\u001b[0;32m--> 161\u001b[0m     derphi1 \u001b[38;5;241m=\u001b[39m \u001b[43mderphi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstp\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    162\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    163\u001b[0m     \u001b[38;5;28;01mbreak\u001b[39;00m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_linesearch.py:78\u001b[0m, in \u001b[0;36mline_search_wolfe1.<locals>.derphi\u001b[0;34m(s)\u001b[0m\n\u001b[1;32m     77\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mderphi\u001b[39m(s):\n\u001b[0;32m---> 78\u001b[0m     gval[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[43mfprime\u001b[49m\u001b[43m(\u001b[49m\u001b[43mxk\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43ms\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mpk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     79\u001b[0m     gc[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m     80\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39mdot(gval[\u001b[38;5;241m0\u001b[39m], pk)\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:273\u001b[0m, in \u001b[0;36mScalarFunction.grad\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    271\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39marray_equal(x, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx):\n\u001b[1;32m    272\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_x_impl(x)\n\u001b[0;32m--> 273\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_update_grad\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    274\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:256\u001b[0m, in \u001b[0;36mScalarFunction._update_grad\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    254\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_update_grad\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m    255\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg_updated:\n\u001b[0;32m--> 256\u001b[0m         \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_update_grad_impl\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    257\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg_updated \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:173\u001b[0m, in \u001b[0;36mScalarFunction.__init__.<locals>.update_grad\u001b[0;34m()\u001b[0m\n\u001b[1;32m    171\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_fun()\n\u001b[1;32m    172\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mngev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m--> 173\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg \u001b[38;5;241m=\u001b[39m \u001b[43mapprox_derivative\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfun_wrapped\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mf0\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    174\u001b[0m \u001b[43m                           \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mfinite_diff_options\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_numdiff.py:505\u001b[0m, in \u001b[0;36mapprox_derivative\u001b[0;34m(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs)\u001b[0m\n\u001b[1;32m    502\u001b[0m     use_one_sided \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[1;32m    504\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m sparsity \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 505\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_dense_difference\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfun_wrapped\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mf0\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mh\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    506\u001b[0m \u001b[43m                             \u001b[49m\u001b[43muse_one_sided\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    507\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    508\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m issparse(sparsity) \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(sparsity) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_numdiff.py:576\u001b[0m, in \u001b[0;36m_dense_difference\u001b[0;34m(fun, x0, f0, h, use_one_sided, method)\u001b[0m\n\u001b[1;32m    574\u001b[0m     x \u001b[38;5;241m=\u001b[39m x0 \u001b[38;5;241m+\u001b[39m h_vecs[i]\n\u001b[1;32m    575\u001b[0m     dx \u001b[38;5;241m=\u001b[39m x[i] \u001b[38;5;241m-\u001b[39m x0[i]  \u001b[38;5;66;03m# Recompute dx as exactly representable number.\u001b[39;00m\n\u001b[0;32m--> 576\u001b[0m     df \u001b[38;5;241m=\u001b[39m \u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m-\u001b[39m f0\n\u001b[1;32m    577\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m method \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m3-point\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m use_one_sided[i]:\n\u001b[1;32m    578\u001b[0m     x1 \u001b[38;5;241m=\u001b[39m x0 \u001b[38;5;241m+\u001b[39m h_vecs[i]\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_numdiff.py:456\u001b[0m, in \u001b[0;36mapprox_derivative.<locals>.fun_wrapped\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    455\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfun_wrapped\u001b[39m(x):\n\u001b[0;32m--> 456\u001b[0m     f \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39matleast_1d(\u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m    457\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m f\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[1;32m    458\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m`fun` return value has \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    459\u001b[0m                            \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmore than 1 dimension.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:137\u001b[0m, in \u001b[0;36mScalarFunction.__init__.<locals>.fun_wrapped\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    133\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnfev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m    134\u001b[0m \u001b[38;5;66;03m# Send a copy because the user may overwrite it.\u001b[39;00m\n\u001b[1;32m    135\u001b[0m \u001b[38;5;66;03m# Overwriting results in undefined behaviour because\u001b[39;00m\n\u001b[1;32m    136\u001b[0m \u001b[38;5;66;03m# fun(self.x) will change self.x, with the two no longer linked.\u001b[39;00m\n\u001b[0;32m--> 137\u001b[0m fx \u001b[38;5;241m=\u001b[39m \u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcopy\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    138\u001b[0m \u001b[38;5;66;03m# Make sure the function returns a true scalar\u001b[39;00m\n\u001b[1;32m    139\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39misscalar(fx):\n",
      "Cell \u001b[0;32mIn[8], line 7\u001b[0m, in \u001b[0;36mEVH\u001b[0;34m(a, psi)\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;66;03m# Calculate normalization of psi\u001b[39;00m\n\u001b[1;32m      6\u001b[0m norm \u001b[38;5;241m=\u001b[39m quad(\u001b[38;5;28;01mlambda\u001b[39;00m x :np\u001b[38;5;241m.\u001b[39mabs(psi(x,a))\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m, xmin, xmax, points\u001b[38;5;241m=\u001b[39mbad_points)[\u001b[38;5;241m0\u001b[39m]\n\u001b[0;32m----> 7\u001b[0m ev \u001b[38;5;241m=\u001b[39m \u001b[43mquad\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43;01mlambda\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m \u001b[49m\u001b[43m:\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mconj\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mHpsi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxmin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mxmax\u001b[49m\u001b[43m,\u001b[49m\u001b[43mpoints\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mbad_points\u001b[49m\u001b[43m)\u001b[49m[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m      8\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ev\u001b[38;5;241m/\u001b[39mnorm\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/integrate/_quadpack_py.py:463\u001b[0m, in \u001b[0;36mquad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst, complex_func)\u001b[0m\n\u001b[1;32m    460\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m retval\n\u001b[1;32m    462\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m weight \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 463\u001b[0m     retval \u001b[38;5;241m=\u001b[39m \u001b[43m_quad\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfull_output\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mepsabs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mepsrel\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlimit\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    464\u001b[0m \u001b[43m                   \u001b[49m\u001b[43mpoints\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    465\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    466\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m points \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/scipy/integrate/_quadpack_py.py:587\u001b[0m, in \u001b[0;36m_quad\u001b[0;34m(func, a, b, args, full_output, epsabs, epsrel, limit, points)\u001b[0m\n\u001b[1;32m    585\u001b[0m the_points \u001b[38;5;241m=\u001b[39m the_points[the_points \u001b[38;5;241m<\u001b[39m b]\n\u001b[1;32m    586\u001b[0m the_points \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mconcatenate((the_points, (\u001b[38;5;241m0.\u001b[39m, \u001b[38;5;241m0.\u001b[39m)))\n\u001b[0;32m--> 587\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_quadpack\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_qagpe\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43mthe_points\u001b[49m\u001b[43m,\u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43mfull_output\u001b[49m\u001b[43m,\u001b[49m\u001b[43mepsabs\u001b[49m\u001b[43m,\u001b[49m\u001b[43mepsrel\u001b[49m\u001b[43m,\u001b[49m\u001b[43mlimit\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[8], line 7\u001b[0m, in \u001b[0;36mEVH.<locals>.<lambda>\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;66;03m# Calculate normalization of psi\u001b[39;00m\n\u001b[1;32m      6\u001b[0m norm \u001b[38;5;241m=\u001b[39m quad(\u001b[38;5;28;01mlambda\u001b[39;00m x :np\u001b[38;5;241m.\u001b[39mabs(psi(x,a))\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m, xmin, xmax, points\u001b[38;5;241m=\u001b[39mbad_points)[\u001b[38;5;241m0\u001b[39m]\n\u001b[0;32m----> 7\u001b[0m ev \u001b[38;5;241m=\u001b[39m quad(\u001b[38;5;28;01mlambda\u001b[39;00m x :np\u001b[38;5;241m.\u001b[39mconj(psi(x,a))\u001b[38;5;241m*\u001b[39m\u001b[43mHpsi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m, xmin, xmax,points\u001b[38;5;241m=\u001b[39mbad_points)[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m      8\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ev\u001b[38;5;241m/\u001b[39mnorm\n",
      "Cell \u001b[0;32mIn[7], line 3\u001b[0m, in \u001b[0;36mHpsi\u001b[0;34m(psi, x, a)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mHpsi\u001b[39m(psi,x,a):\n\u001b[1;32m      2\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m'''The Hamiltonian acting on psi at point x'''\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m0.5\u001b[39m\u001b[38;5;241m*\u001b[39mhbar\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m/\u001b[39mmass \u001b[38;5;241m*\u001b[39m\u001b[43mnd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mDerivative\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43;01mlambda\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m+\u001b[39m V(x)\u001b[38;5;241m*\u001b[39mpsi(x,a)\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numdifftools/core.py:288\u001b[0m, in \u001b[0;36mDerivative.__call__\u001b[0;34m(self, x, *args, **kwds)\u001b[0m\n\u001b[1;32m    286\u001b[0m x_i \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(x)\n\u001b[1;32m    287\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m np\u001b[38;5;241m.\u001b[39merrstate(divide\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mignore\u001b[39m\u001b[38;5;124m'\u001b[39m, invalid\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mignore\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m--> 288\u001b[0m     results, f_xi \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_derivative\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx_i\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    289\u001b[0m     derivative, info \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_extrapolate(\u001b[38;5;241m*\u001b[39mresults)\n\u001b[1;32m    290\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfull_output:\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numdifftools/core.py:238\u001b[0m, in \u001b[0;36mDerivative._derivative_nonzero_order\u001b[0;34m(self, x_i, args, kwds)\u001b[0m\n\u001b[1;32m    236\u001b[0m steps, step_ratio \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_steps(x_i)\n\u001b[1;32m    237\u001b[0m fxi \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_eval_first(f, x_i)\n\u001b[0;32m--> 238\u001b[0m results \u001b[38;5;241m=\u001b[39m [diff(f, fxi, x_i, h) \u001b[38;5;28;01mfor\u001b[39;00m h \u001b[38;5;129;01min\u001b[39;00m steps]\n\u001b[1;32m    240\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mset_richardson_rule(step_ratio, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrichardson_terms)\n\u001b[1;32m    242\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfd_rule\u001b[38;5;241m.\u001b[39mapply(results, steps, step_ratio), fxi\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numdifftools/core.py:238\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    236\u001b[0m steps, step_ratio \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_steps(x_i)\n\u001b[1;32m    237\u001b[0m fxi \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_eval_first(f, x_i)\n\u001b[0;32m--> 238\u001b[0m results \u001b[38;5;241m=\u001b[39m [\u001b[43mdiff\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfxi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx_i\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mh\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m h \u001b[38;5;129;01min\u001b[39;00m steps]\n\u001b[1;32m    240\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mset_richardson_rule(step_ratio, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrichardson_terms)\n\u001b[1;32m    242\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfd_rule\u001b[38;5;241m.\u001b[39mapply(results, steps, step_ratio), fxi\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numdifftools/finite_difference.py:57\u001b[0m, in \u001b[0;36mDifferenceFunctions._central_even\u001b[0;34m(f, f_x0i, x0i, h)\u001b[0m\n\u001b[1;32m     55\u001b[0m \u001b[38;5;129m@staticmethod\u001b[39m\n\u001b[1;32m     56\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_central_even\u001b[39m(f, f_x0i, x0i, h):\n\u001b[0;32m---> 57\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m (\u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx0i\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mh\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m+\u001b[39m f(x0i \u001b[38;5;241m-\u001b[39m h)) \u001b[38;5;241m/\u001b[39m \u001b[38;5;241m2.0\u001b[39m \u001b[38;5;241m-\u001b[39m f_x0i\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numdifftools/core.py:257\u001b[0m, in \u001b[0;36mDerivative._get_functions.<locals>.export_fun\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    256\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mexport_fun\u001b[39m(x):\n\u001b[0;32m--> 257\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfun\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[7], line 3\u001b[0m, in \u001b[0;36mHpsi.<locals>.<lambda>\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mHpsi\u001b[39m(psi,x,a):\n\u001b[1;32m      2\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m'''The Hamiltonian acting on psi at point x'''\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m0.5\u001b[39m\u001b[38;5;241m*\u001b[39mhbar\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m/\u001b[39mmass \u001b[38;5;241m*\u001b[39mnd\u001b[38;5;241m.\u001b[39mDerivative(\u001b[38;5;28;01mlambda\u001b[39;00m x: \u001b[43mpsi\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m, n\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m)(x) \u001b[38;5;241m+\u001b[39m V(x)\u001b[38;5;241m*\u001b[39mpsi(x,a)\n",
      "Cell \u001b[0;32mIn[5], line 7\u001b[0m, in \u001b[0;36mpsi\u001b[0;34m(x, a)\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m'''Wavefunction as a function of position x.\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;124;03mThe array \"a\" stores any parameters.'''\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;66;03m# Note, normalization is not necessary\u001b[39;00m\n\u001b[1;32m      6\u001b[0m \u001b[38;5;66;03m#return np.exp(-x**2*a[0]) # gaussian\u001b[39;00m\n\u001b[0;32m----> 7\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpiecewise\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m[\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mabs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43ma\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m<\u001b[39;49m\u001b[38;5;241;43m1.0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mabs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43ma\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m>\u001b[39;49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1.0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m1.0\u001b[39;49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcos\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpi\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m0.0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36mpiecewise\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numpy/lib/function_base.py:749\u001b[0m, in \u001b[0;36mpiecewise\u001b[0;34m(x, condlist, funclist, *args, **kw)\u001b[0m\n\u001b[1;32m    743\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m n \u001b[38;5;241m!=\u001b[39m n2:\n\u001b[1;32m    744\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    745\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mwith \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m condition(s), either \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m or \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m functions are expected\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    746\u001b[0m         \u001b[38;5;241m.\u001b[39mformat(n, n, n\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m    747\u001b[0m     )\n\u001b[0;32m--> 749\u001b[0m y \u001b[38;5;241m=\u001b[39m \u001b[43mzeros_like\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    750\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m cond, func \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(condlist, funclist):\n\u001b[1;32m    751\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(func, collections\u001b[38;5;241m.\u001b[39mabc\u001b[38;5;241m.\u001b[39mCallable):\n",
      "File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36mzeros_like\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numpy/core/numeric.py:142\u001b[0m, in \u001b[0;36mzeros_like\u001b[0;34m(a, dtype, order, subok, shape)\u001b[0m\n\u001b[1;32m    140\u001b[0m \u001b[38;5;66;03m# needed instead of a 0 to get same result as zeros for string dtypes\u001b[39;00m\n\u001b[1;32m    141\u001b[0m z \u001b[38;5;241m=\u001b[39m zeros(\u001b[38;5;241m1\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mres\u001b[38;5;241m.\u001b[39mdtype)\n\u001b[0;32m--> 142\u001b[0m \u001b[43mmultiarray\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcopyto\u001b[49m\u001b[43m(\u001b[49m\u001b[43mres\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mz\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcasting\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43munsafe\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    143\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m res\n",
      "File \u001b[0;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36mcopyto\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[0;32m~/mambaforge/envs/scipy/lib/python3.10/site-packages/numpy/core/multiarray.py:1079\u001b[0m, in \u001b[0;36mcopyto\u001b[0;34m(dst, src, casting, where)\u001b[0m\n\u001b[1;32m   1031\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1032\u001b[0m \u001b[38;5;124;03m    unravel_index(indices, shape, order='C')\u001b[39;00m\n\u001b[1;32m   1033\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1074\u001b[0m \n\u001b[1;32m   1075\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[1;32m   1076\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m (indices,)\n\u001b[0;32m-> 1079\u001b[0m \u001b[38;5;129m@array_function_from_c_func_and_dispatcher\u001b[39m(_multiarray_umath\u001b[38;5;241m.\u001b[39mcopyto)\n\u001b[1;32m   1080\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcopyto\u001b[39m(dst, src, casting\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, where\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m   1081\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1082\u001b[0m \u001b[38;5;124;03m    copyto(dst, src, casting='same_kind', where=True)\u001b[39;00m\n\u001b[1;32m   1083\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   1125\u001b[0m \u001b[38;5;124;03m       \u001b[39;00m\n\u001b[1;32m   1126\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[1;32m   1127\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m (dst, src, where)\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# minimize the expectation value of energy by changing a_guess\n",
    "# `output` has two pieces\n",
    "# `output.x` is the list of parameters (what I called a above)\n",
    "output = minimize(EVH, a_guess, args=(psi,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:02.723920Z",
     "start_time": "2018-10-19T18:44:02.708301Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# build arrays for plotting x, V(x), and psi(x)\n",
    "# these are only used for plotting\n",
    "xlist = np.linspace(xmin, xmax, 101)\n",
    "ylist = np.array([V(x) for x in xlist])\n",
    "\n",
    "# intital and final wavefunctions\n",
    "psilist = np.array([psi(x,a_guess) for x in xlist])\n",
    "outlist = np.array([psi(x,output.x) for x in xlist])\n",
    "\n",
    "# probability density functions of t\n",
    "problist = np.array([np.abs(psi(x,output.x))**2 for x in xlist])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:03.192655Z",
     "start_time": "2018-10-19T18:44:03.177033Z"
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(\"Minimum energy: {} for arguments {}\".format(output.fun,output.x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-19T18:44:07.145659Z",
     "start_time": "2018-10-19T18:44:06.864416Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(xlist,ylist,'k-', label = \"V(x)\")\n",
    "plt.plot(xlist,(psilist+1)*initEV,'b-', label = \"Initial ψ(x)\")\n",
    "plt.plot(xlist,(outlist+1)*output.fun,'r-', label = \"Final ψ(x)\")\n",
    "#plt.plot(xlist,problist+output.fun, 'g-', label = \"|ψ(x)|^2\")\n",
    "plt.legend();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"Energy or ψ\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian psi: Minimum energy: 0.6970861111798174 for arguments 1.87843196"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (scipy)",
   "language": "python",
   "name": "scipy"
  },
  "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.10.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}