{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "cells": [
  {
   "cell_type": "code",
   "id": "cell-0",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# --- house style (the Press palette: accent/navy/gold/parchment) ---\nimport matplotlib as mpl\nACCENT, NAVY, GOLD, PARCH = '#7a1f1f', '#1f3a5f', '#b8860b', '#f7f2e7'\nmpl.rcParams.update({'figure.facecolor': PARCH, 'axes.facecolor': '#fffdf6',\n                     'axes.edgecolor': '#c9bfa3', 'font.family': 'serif'})\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nRULES = {'bump':       lambda x: np.exp(-((x-0.5)**2)/0.05) + 0.25,\n         'wall':       lambda x: 1.0/(0.05 + x),\n         'two peaks':  lambda x: np.exp(-((x-0.3)**2)/0.01) + 0.7*np.exp(-((x-0.75)**2)/0.005) + 0.1,\n         'constant':   lambda x: np.full_like(x, 2.0)}     # the equitesseral check\nrho = RULES['bump']                            # <- choose your tiling rule here\n\nxx = np.linspace(0, 1, 2001)\nmu = np.trapezoid(rho(xx), xx)\nV = np.array([np.trapezoid(rho(xx[:i+1]), xx[:i+1]) for i in range(len(xx))]) / mu\ndV = np.gradient(V, xx)\ntheory = rho(xx) / mu\n\nx0 = 0.62\ni0 = np.argmin(np.abs(xx - x0))\nprint(f\"at x = {x0}:  numeric dV/dx = {dV[i0]:.5f}   rho(x)/mu = {theory[i0]:.5f}  (Thm 10.2)\")\nprint(f\"worst |numeric - theory| away from edges: \"\n      f\"{np.abs(dV - theory)[50:-50].max():.2e}\")\n\nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)\nax1.plot(xx, rho(xx), color=NAVY, lw=2)\nax1.fill_between(xx[xx <= x0], rho(xx[xx <= x0]), color=ACCENT, alpha=0.4)\nax1.set(ylabel='rho(x)', title='the tiling rule, swept up to x')\nax2.plot(xx, V, color=ACCENT, lw=2, label='V(x): the running tesseraction')\nax2.plot(xx, dV, color=GOLD, lw=1.5, ls='--', label=\"numeric V'(x)\")\nax2.plot(xx, theory, color=NAVY, lw=1.2, ls=':', label='rho(x)/mu')\nax2.plot(x0, V[i0], 'o', color='k')\nax2.set(xlabel='x', ylabel='share', title='differentiation and integration cancel')\nax2.legend()\nplt.tight_layout(); plt.show()\n"
  }
 ]
}