Lecture 13.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations\n",
    "gr();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# M3M6: Methods of Mathematical Physics\n",
    "\n",
    "$$\n",
    "\\def\\dashint{{\\int\\!\\!\\!\\!\\!\\!-\\,}}\n",
    "\\def\\infdashint{\\dashint_{\\!\\!\\!-\\infty}^{\\,\\infty}}\n",
    "\\def\\D{\\,{\\rm d}}\n",
    "\\def\\dx{\\D x}\n",
    "\\def\\dt{\\D t}\n",
    "\\def\\C{{\\mathbb C}}\n",
    "\\def\\CC{{\\cal C}}\n",
    "\\def\\HH{{\\cal H}}\n",
    "\\def\\I{{\\rm i}}\n",
    "\\def\\qqfor{\\qquad\\hbox{for}\\qquad}\n",
    "$$\n",
    "\n",
    "Dr. Sheehan Olver\n",
    "<br>\n",
    "s.olver@imperial.ac.uk\n",
    "\n",
    "<br>\n",
    "Website: https://github.com/dlfivefifty/M3M6LectureNotes\n",
    "\n",
    "# Lecture 13: Electric charges in a potential well\n",
    "\n",
    "1. One charge in a well\n",
    "2. Two charges in a well\n",
    "3. $N$ charges in a well\n",
    "4. Limiting distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We study the dynamics of many electric charges in a potential well.  We restrict our attention to 1D: picture an infinitely long wire with carges on it. We will see that as the number of charges becomes large, we can determine the limiting distribution using an additive Riemann–Hilbert problem, but one which the _interval_ its posed on is solved for.\n",
    "\n",
    "## One charge in a potential well \n",
    "\n",
    "Consider a point charge in a well $V(x) = x^2 / 2$, initially located at $\\lambda^0$. The dynamics of the point charge are governed by\n",
    "$$\n",
    "{ \\D \\lambda \\over \\dt} = -V'(\\lambda) = -  \\lambda\n",
    "$$\n",
    "that is: if we are positive we move left and if we are negative we move right. Here is a movie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = x -> x^2/2\n",
    "Vp = x -> x\n",
    "\n",
    "λ⁰ = 2.3  # initial location\n",
    "prob = ODEProblem((λ,_,t) -> -Vp(λ), λ⁰, (0.0, 10.0))\n",
    "λ = solve(prob; reltol=1E-6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "@gif for t=0.0:0.05:7.0\n",
    "    plot(V, range(-5, stop=5, length=100); label=\"potential\", title=\"t = $t\")\n",
    "    scatter!([λ(t)] ,[0.0]; label=\"charge\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![One charge](onecharge.gif)\n",
    "\n",
    "In the limit the charge reaches an equilibrium: it no longer varies in time. I.e., it reaches a point where ${\\D\\lambda \\over \\dt} = 0$, which is equivalent to solving\n",
    "$$\n",
    " 0 = - V'(\\lambda)  = - \\lambda\n",
    "$$\n",
    "in other words, the minimum of the well, in this case $\\lambda = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two charges in a potential well \n",
    "\n",
    "Suppose there are now two charges, $\\lambda_1$ and $\\lambda_2$.  The effect on the first charge $\\lambda_1$ is to repulse away from $\\lambda_2$ via via:\n",
    "$$\n",
    "{\\D \\lambda_1 \\over \\D t} = {1 \\over \\lambda_1 -\\lambda_2}\n",
    "$$\n",
    "Similarly, the effect on $\\lambda_2$ is\n",
    "$$\n",
    "{\\D \\lambda_2 \\over \\D t} = {1 \\over \\lambda_2 -\\lambda_1}\n",
    "$$\n",
    "Unrestricted, the two potentials will repulse off to infinity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ⁰ = [1.2,2.3]  # initial location\n",
    "prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]), 1/(λ[2] - λ[1])], λ⁰, (0.0, 10.0))\n",
    "λ = solve(prob; reltol=1E-6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@gif for t=0.0:0.05:7.0\n",
    "    scatter(λ(t) ,zeros(2); label=\"charges\", xlims=(-5,5), title=\"t = $t\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Two charges, no potential](twocharges_nov.gif)\n",
    "\n",
    "Adding in a potential well and we get an equilbrium again:\n",
    "\\begin{align*}\n",
    "{\\D \\lambda_1 \\over \\D t} = {1 \\over \\lambda_1 -\\lambda_2} - V'(\\lambda_1) \\\\\n",
    "{\\D \\lambda_2 \\over \\D t} = {1 \\over \\lambda_2 -\\lambda_1} - V'(\\lambda_2)\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ₀ = [1.2,2.3]  # initial location\n",
    "prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]) - Vp(λ[1]), 1/(λ[2] - λ[1]) - Vp(λ[2])], λ₀, (0.0, 10.0))\n",
    "λ = solve(prob; reltol=1E-6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@gif for t=0.0:0.05:7.0\n",
    "    plot(V, range(-5, stop=5, length=100); label=\"potential\", title=\"t = $t\")\n",
    "    scatter!(λ(t) ,zeros(2); label=\"charges\", xlims=(-5,5), title=\"t = $t\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Two charges](twocharges.gif)\n",
    "\n",
    "The limiting distribution is given by\n",
    "\\begin{align*}\n",
    "0 = {1 \\over \\lambda_1 -\\lambda_2} - V'(\\lambda_1) \\\\\n",
    "0 = {1 \\over \\lambda_2 -\\lambda_1} - V'(\\lambda_2) \n",
    "\\end{align*}\n",
    "For this potential, we can solve it exactly: we need to solve\n",
    "\\begin{align*}\n",
    "\\lambda_1 = {1 \\over \\lambda_1 -\\lambda_2}  \\\\\n",
    "\\lambda_2 = {1 \\over \\lambda_2 -\\lambda_1} \n",
    "\\end{align*}\n",
    "Using $\\lambda_1 = -\\lambda_2$, we find that $\\lambda_1 = \\pm{1 \\over \\sqrt 2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $N$ charges in a potential well\n",
    "\n",
    "Each charge repulses every other charge, so we end up needing to sum over them all:\n",
    "\\begin{align*}\n",
    "{\\D \\lambda_k \\over \\D t} = \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\lambda_k -\\lambda_j} - V'(\\lambda_k) \\\\\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "λ⁰ = randn(N)  # initial location\n",
    "prob = ODEProblem((λ,_,t) -> [sum(1 ./ (λ[k] .- λ[[1:k-1;k+1:end]])) - Vp(λ[k]) for k=1:N], λ⁰, (0.0, 10.0))\n",
    "λ = solve(prob; reltol=1E-6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@gif for t=0.0:0.05:7.0\n",
    "    plot(V, range(-5, stop=5, length=100); label=\"potential\", title=\"t = $t\")\n",
    "    scatter!(λ(t) ,zeros(N); label=\"charges\", xlims=(-5,5), title=\"t = $t\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![N charges](ncharges.gif)\n",
    "\n",
    "As the number of charges becomes large, they spread off to infinity. In the case of $V(x) = x^2$, we can renormalize by  dividing by $N$ so they stay bounded:\n",
    "$\\mu_k = {\\lambda_k \\over \\sqrt N}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@gif for t=0.0:0.05:7.0\n",
    "    scatter(λ(t)/sqrt(N) ,zeros(N); label=\"charges\", xlims=(-2,2), title=\"t = $t\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![N charges, rescaled](ncharges_scaled.gif)\n",
    "\n",
    "This begs questions: why does it balance out at $\\pm \\sqrt 2$? Why does it have a nice histogram precisely like ${\\sqrt{2-x^2} \\over \\pi}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip1400\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip1401\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip1401)\" points=\"\n",
       "0,1600 2400,1600 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip1402\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip1401)\" points=\"\n",
       "161.394,1503.47 2321.26,1503.47 2321.26,47.2441 161.394,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip1403\">\n",
       "    <rect x=\"161\" y=\"47\" width=\"2161\" height=\"1457\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  211.538,1503.47 211.538,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  554.801,1503.47 554.801,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  898.064,1503.47 898.064,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1241.33,1503.47 1241.33,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1584.59,1503.47 1584.59,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1927.85,1503.47 1927.85,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2271.12,1503.47 2271.12,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1462.26 2321.26,1462.26 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1157.06 2321.26,1157.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,851.865 2321.26,851.865 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,546.667 2321.26,546.667 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,241.47 2321.26,241.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1503.47 2321.26,1503.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1503.47 161.394,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  211.538,1503.47 211.538,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  554.801,1503.47 554.801,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  898.064,1503.47 898.064,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1241.33,1503.47 1241.33,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1584.59,1503.47 1584.59,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1927.85,1503.47 1927.85,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2271.12,1503.47 2271.12,1481.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1462.26 193.792,1462.26 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1157.06 193.792,1157.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,851.865 193.792,851.865 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,546.667 193.792,546.667 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,241.47 193.792,241.47 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 211.538, 1557.47)\" x=\"211.538\" y=\"1557.47\">-1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 554.801, 1557.47)\" x=\"554.801\" y=\"1557.47\">-1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 898.064, 1557.47)\" x=\"898.064\" y=\"1557.47\">-0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1241.33, 1557.47)\" x=\"1241.33\" y=\"1557.47\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1584.59, 1557.47)\" x=\"1584.59\" y=\"1557.47\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1927.85, 1557.47)\" x=\"1927.85\" y=\"1557.47\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2271.12, 1557.47)\" x=\"2271.12\" y=\"1557.47\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1479.76)\" x=\"137.394\" y=\"1479.76\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1174.56)\" x=\"137.394\" y=\"1174.56\">0.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 869.365)\" x=\"137.394\" y=\"869.365\">0.2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 564.167)\" x=\"137.394\" y=\"564.167\">0.3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 258.97)\" x=\"137.394\" y=\"258.97\">0.4</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "280.191,1004.46 280.191,1462.26 417.496,1462.26 417.496,1004.46 280.191,1004.46 280.191,1004.46 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  280.191,1004.46 280.191,1462.26 417.496,1462.26 417.496,1004.46 280.191,1004.46 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "417.496,546.667 417.496,1462.26 554.801,1462.26 554.801,546.667 417.496,546.667 417.496,546.667 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  417.496,546.667 417.496,1462.26 554.801,1462.26 554.801,546.667 417.496,546.667 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "554.801,394.069 554.801,1462.26 692.106,1462.26 692.106,394.069 554.801,394.069 554.801,394.069 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  554.801,394.069 554.801,1462.26 692.106,1462.26 692.106,394.069 554.801,394.069 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "692.106,241.47 692.106,1462.26 829.411,1462.26 829.411,241.47 692.106,241.47 692.106,241.47 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  692.106,241.47 692.106,1462.26 829.411,1462.26 829.411,241.47 692.106,241.47 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "829.411,241.47 829.411,1462.26 966.717,1462.26 966.717,241.47 829.411,241.47 829.411,241.47 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  829.411,241.47 829.411,1462.26 966.717,1462.26 966.717,241.47 829.411,241.47 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "966.717,88.8708 966.717,1462.26 1104.02,1462.26 1104.02,88.8708 966.717,88.8708 966.717,88.8708 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  966.717,88.8708 966.717,1462.26 1104.02,1462.26 1104.02,88.8708 966.717,88.8708 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1104.02,88.8708 1104.02,1462.26 1241.33,1462.26 1241.33,88.8708 1104.02,88.8708 1104.02,88.8708 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1104.02,88.8708 1104.02,1462.26 1241.33,1462.26 1241.33,88.8708 1104.02,88.8708 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1241.33,88.8708 1241.33,1462.26 1378.63,1462.26 1378.63,88.8708 1241.33,88.8708 1241.33,88.8708 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1241.33,88.8708 1241.33,1462.26 1378.63,1462.26 1378.63,88.8708 1241.33,88.8708 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1378.63,88.8708 1378.63,1462.26 1515.94,1462.26 1515.94,88.8708 1378.63,88.8708 1378.63,88.8708 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1378.63,88.8708 1378.63,1462.26 1515.94,1462.26 1515.94,88.8708 1378.63,88.8708 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1515.94,241.47 1515.94,1462.26 1653.24,1462.26 1653.24,241.47 1515.94,241.47 1515.94,241.47 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1515.94,241.47 1515.94,1462.26 1653.24,1462.26 1653.24,241.47 1515.94,241.47 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1653.24,241.47 1653.24,1462.26 1790.55,1462.26 1790.55,241.47 1653.24,241.47 1653.24,241.47 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1653.24,241.47 1653.24,1462.26 1790.55,1462.26 1790.55,241.47 1653.24,241.47 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1790.55,394.069 1790.55,1462.26 1927.85,1462.26 1927.85,394.069 1790.55,394.069 1790.55,394.069 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1790.55,394.069 1790.55,1462.26 1927.85,1462.26 1927.85,394.069 1790.55,394.069 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "1927.85,546.667 1927.85,1462.26 2065.16,1462.26 2065.16,546.667 1927.85,546.667 1927.85,546.667 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1927.85,546.667 1927.85,1462.26 2065.16,1462.26 2065.16,546.667 1927.85,546.667 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1403)\" points=\"\n",
       "2065.16,1004.46 2065.16,1462.26 2202.46,1462.26 2202.46,1004.46 2065.16,1004.46 2065.16,1004.46 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2065.16,1004.46 2065.16,1462.26 2202.46,1462.26 2202.46,1004.46 2065.16,1004.46 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1403)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.433,1462.26 290.047,1187.5 309.661,1075.68 329.275,991.243 348.889,921.216 368.503,860.547 388.117,806.631 407.731,757.918 427.345,713.39 446.959,672.339 \n",
       "  466.573,634.249 486.187,598.728 505.801,565.469 525.415,534.23 545.029,504.814 564.643,477.056 584.257,450.82 603.871,425.991 623.485,402.47 643.099,380.171 \n",
       "  662.713,359.021 682.327,338.954 701.941,319.914 721.555,301.849 741.169,284.716 760.783,268.473 780.397,253.085 800.011,238.52 819.625,224.749 839.239,211.745 \n",
       "  858.853,199.484 878.467,187.946 898.081,177.111 917.695,166.961 937.309,157.48 956.924,148.654 976.538,140.47 996.152,132.915 1015.77,125.979 1035.38,119.653 \n",
       "  1054.99,113.927 1074.61,108.795 1094.22,104.25 1113.84,100.285 1133.45,96.8951 1153.06,94.077 1172.68,91.8267 1192.29,90.1414 1211.91,89.019 1231.52,88.4582 \n",
       "  1251.13,88.4582 1270.75,89.019 1290.36,90.1414 1309.98,91.8267 1329.59,94.077 1349.2,96.8951 1368.82,100.285 1388.43,104.25 1408.05,108.795 1427.66,113.927 \n",
       "  1447.27,119.653 1466.89,125.979 1486.5,132.915 1506.12,140.47 1525.73,148.654 1545.34,157.48 1564.96,166.961 1584.57,177.111 1604.19,187.946 1623.8,199.484 \n",
       "  1643.41,211.745 1663.03,224.749 1682.64,238.52 1702.26,253.085 1721.87,268.473 1741.48,284.716 1761.1,301.849 1780.71,319.914 1800.33,338.954 1819.94,359.021 \n",
       "  1839.55,380.171 1859.17,402.47 1878.78,425.991 1898.4,450.82 1918.01,477.056 1937.62,504.814 1957.24,534.23 1976.85,565.469 1996.47,598.728 2016.08,634.249 \n",
       "  2035.7,672.339 2055.31,713.39 2074.92,757.918 2094.54,806.631 2114.15,860.547 2133.77,921.216 2153.38,991.243 2172.99,1075.68 2192.61,1187.5 2212.22,1462.26 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1401)\" points=\"\n",
       "1559.84,312.204 2249.26,312.204 2249.26,130.764 1559.84,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1559.84,312.204 2249.26,312.204 2249.26,130.764 1559.84,130.764 1559.84,312.204 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip1401)\" points=\"\n",
       "1583.84,215.436 1727.84,215.436 1727.84,167.052 1583.84,167.052 1583.84,215.436 \n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1583.84,215.436 1727.84,215.436 1727.84,167.052 1583.84,167.052 1583.84,215.436 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1751.84, 208.744)\" x=\"1751.84\" y=\"208.744\">histogram of charges</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip1401)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1583.84,251.724 1727.84,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip1401)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1751.84, 269.224)\" x=\"1751.84\" y=\"269.224\">semicircle</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "histogram(λ(10.0)/sqrt(N); nbins=20, normalize=true, label=\"histogram of charges\")\n",
    "plot!(x -> sqrt(2-x^2)/(π), range(eps()-sqrt(2.0); stop=sqrt(2)-eps(), length=100), label=\"semicircle\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equilibrium distribution\n",
    "\n",
    "Plugging in $\\lambda_k = \\sqrt N \\mu_k$, we get a dynamical system for $\\mu_k$:\n",
    "$$\n",
    " {\\D  \\mu_k \\over \\D t} = {1 \\over N} \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\mu_k -\\mu_j} -  \\mu_k\n",
    "$$\n",
    "(The choice of scaling like $\\sqrt N$ was dictated by $V'(x)$, if $V(x) = x^4$ it would be $N^{1/4}$.) Thus the limit of the charges is given by\n",
    "$$\n",
    "0 = {1 \\over N} \\sum_{j=1 \\atop j \\neq k}^N {1 \\over \\mu_k -\\mu_j} -  \\mu_k \n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "It is convenient to represent the point charges by Dirac delta functions:\n",
    "$$\n",
    "    w_N(x) = {1 \\over N} \\sum_{k=1}^N \\delta_{\\mu_k}(x)\n",
    "$$\n",
    "normalized so that $\\int w_N(x) \\dx = 1$, so that\n",
    "$$\n",
    "    {1 \\over N} \\sum_{k=1}^N {1 \\over  x -\\mu_j} = \\int_{-\\infty}^\\infty {w_N(t) \\dt \\over x - t}\n",
    "$$\n",
    "or in other words, we have \n",
    "$$\n",
    " \\HH_{(-\\infty,\\infty)} w_N(\\mu_k) = -{V'(\\mu_k) \\over \\pi}\n",
    "$$\n",
    "since\n",
    "$$\n",
    "\\HH w_n (\\mu_k) = {1 \\over \\pi} \\lim_{\\epsilon\\rightarrow 0} \\left(\\int_{-\\infty}^{\\mu_k-\\epsilon} + \\int_{\\mu_k+\\epsilon}^\\infty\\right) {w_N(t)  \\over t - \\mu_k} \\dt  = {1 \\over N\\pi} \\sum_{j \\neq k} {1 \\over \\mu_j - \\mu_k}\n",
    "$$\n",
    "\n",
    "Formally (see a more detailed explanation below), $w_N(x)$ tends to a continuous limit as $N\\rightarrow \\infty$, which we have guessed from the histogram to be $w(x) = { \\sqrt{2-x^2}  \\over \\pi}$ for $-\\sqrt 2 < x < \\sqrt2$.  We expect this limit to satisfy the same equation as $w_n$, that is\n",
    "$$\n",
    " \\HH w(x) = -{x \\over \\pi}\n",
    "$$\n",
    "for $x$ in the support of $w(x)$.\n",
    "\n",
    "\n",
    "Indeed:\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: norm not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: norm not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[16]:3"
     ]
    }
   ],
   "source": [
    "x = Fun(-sqrt(2) .. sqrt(2))\n",
    "w = sqrt(2-x^2)/π\n",
    "\n",
    "norm(hilbert(w) + x/π)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is it $[-\\sqrt 2, \\sqrt 2]$? \n",
    "\n",
    "We thus want to choose the interval $[a,b]$ so that there exists a $w(x)$  satisfying\n",
    "\n",
    "1. $w$ is bounded (Based on observation)\n",
    "2. $w(x) \\geq 0$ for $a \\leq x \\leq b$ (Since it is a probability distribution)\n",
    "2. $\\int_a^b w(x) \\dx = 1$ (Since it is a probability distribution)\n",
    "3. $\\HH_{[a,b]} w(x) = - x/\\pi$ \n",
    "\n",
    "As we saw last lecture, there exists a bounded solution to \n",
    "$\\HH_{[-b,b]} u =  - x/\\pi$,\n",
    "namely  $u(x) = { \\sqrt{b^2-x^2} \\over \\pi}$. The choice $b = \\sqrt{2}$ ensures that $\\int_{-b}^b u(x) \\dx = 1$, hence $u(x) = w(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside: Explanation of limit of $w_N(x)$\n",
    "\n",
    "This is beyond the scope of the course, but the convergence of $w_N(x)$ to $w(x)$ is known as weak-\\* convergence. A simple version of this is that\n",
    "$$\n",
    "   \\int_c^d w_N(x) \\dx \\rightarrow \\int_c^d w(x) \\dx\n",
    "$$\n",
    "for every choice of interval $(c,d)$. $   \\int_c^d w_N(x) \\dx $ is precisely the number of charges in $(c,d)$ scaled by $1/N$,  which is exactly what a histogram plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1790059168895313"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = -0.1; b= 0.3;\n",
    "w = Fun(x -> sqrt(2-x^2)/π, a .. b)\n",
    "sum(w)  # integral of w(x) between a and b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.18"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length(filter(λ -> a ≤ λ ≤ b, λ(10.0)/sqrt(N)))/N # integral of w_N(x) between a and b "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another varient of describing weak-\\* convergence is that the Cauchy transforms converge, that is, for $z$ on any  contour $\\gamma$ surrounding $a$ and $b$ (now the support of $w$) we have\n",
    "$$\n",
    "\\int_a^b {w_N(x) \\over x - z} \\dx \\rightarrow \\int_a^b {w(x) \\over x - z} \\dx\n",
    "$$\n",
    "converges uniformly with respect to $N \\rightarrow \\infty$.\n",
    "here we demonstrate it converges at a point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1949409042727309 + 0.013681505291622225im"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Fun(-sqrt(2) .. sqrt(2))\n",
    "w = sqrt(2-x^2)/π\n",
    "z = 0.1+0.2im\n",
    "cauchy(w, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.19542123392652283 + 0.013724448552242337im"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(sum(1 ./((λ(10.0)/sqrt(N)) .- z))/N)/(2π*im)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.2",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}