{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a22f0542",
   "metadata": {},
   "source": [
    "# Laplacian in Spherical Coordinates\n",
    "\n",
    "The Euclidean Laplacian has a coordinate-independent meaning, but\n",
    "its formula changes with the coordinates.  We derive the\n",
    "spherical-coordinate expression from the metric and then\n",
    "verify it by the multivariable chain rule.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0149855e",
   "metadata": {},
   "source": [
    "## Coordinate metric\n",
    "\n",
    "The line element is\n",
    "\n",
    "$$ds^2=dr^2+r^2d\\theta^2+r^2\\sin^2\\theta\\,d\\phi^2.$$\n",
    "\n",
    "The coordinate map is differentiated to obtain tangent vectors. Their dot products give the metric, and its inverse is then computed from that induced matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb5e2dda",
   "metadata": {},
   "outputs": [],
   "source": [
    "declare symbol r, θ, φ : MathValue\n",
    "\n",
    "def x : Vector MathValue := [| r, θ, φ |]\n",
    "def position : Vector MathValue :=\n",
    "  [| r * sin θ * cos φ\n",
    "   , r * sin θ * sin φ\n",
    "   , r * cos θ\n",
    "   |]\n",
    "\n",
    "def e_i_j : Matrix MathValue := ∂/∂ position_j x~i\n",
    "def g_i_j : Matrix MathValue :=\n",
    "  generateTensor (\\[a, b] -> V.* e_a_# e_b_#) [3, 3]\n",
    "def g~i~j : Matrix MathValue := M.inverse g_#_#\n",
    "\n",
    "def f : MathValue := function (r, θ, φ)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bc1e101a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "$\\begin{pmatrix} 1 & 0 & 0 \\\\ 0 & r^{2} & 0 \\\\ 0 & 0 & \\sin(θ)^{2} r^{2} \\\\ \\end{pmatrix}_{\\#\\#}^{\\;\\;}$"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g_#_#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "407069d9",
   "metadata": {},
   "source": [
    "## Covariant derivation\n",
    "\n",
    "For a scalar field, the Laplace--Beltrami operator can be written\n",
    "\n",
    "$$\n",
    "\\Delta f=g^{ij}\\partial_i\\partial_jf\n",
    "   -g^{ij}\\Gamma^k{}_{ij}\\partial_kf.\n",
    "$$\n",
    "\n",
    "Egison contracts the repeated tensor indices in precisely this\n",
    "formula.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "56f9946b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Γ_i_j_k : Tensor MathValue :=\n",
    "  (1 / 2) *\n",
    "    (∂/∂ g_i_k x~j + ∂/∂ g_i_j x~k - ∂/∂ g_j_k x~i)\n",
    "\n",
    "def Γ~i_j_k : Tensor MathValue := withSymbols [m]\n",
    "  g~i~m . Γ_m_j_k\n",
    "\n",
    "def laplacian : MathValue := withSymbols [i, j, k]\n",
    "  g~i~j . ∂/∂ (∂/∂ f x~j) x~i\n",
    "    - g~i~j . Γ~k_i_j . ∂/∂ f x~k\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b628bf7c",
   "metadata": {},
   "source": [
    "## Result\n",
    "\n",
    "$$\n",
    "\\Delta f=\\frac{\\partial^2f}{\\partial r^2}\n",
    "  +\\frac2r\\frac{\\partial f}{\\partial r}\n",
    "  +\\frac1{r^2}\\frac{\\partial^2f}{\\partial\\theta^2}\n",
    "  +\\frac{\\cos\\theta}{r^2\\sin\\theta}\\frac{\\partial f}{\\partial\\theta}\n",
    "  +\\frac1{r^2\\sin^2\\theta}\\frac{\\partial^2f}{\\partial\\phi^2}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2a4bc076",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "$\\frac{\\partial^2 f}{\\partial 1^2} + 2 \\frac{\\partial f}{\\partial 1} r^{-1} + \\frac{\\partial^2 f}{\\partial 2^2} r^{-2} + \\cos(θ) \\sin(θ)^{-1} \\frac{\\partial f}{\\partial 2} r^{-2} + \\sin(θ)^{-2} \\frac{\\partial^2 f}{\\partial 3^2} r^{-2}$"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "laplacian\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81e81f98",
   "metadata": {},
   "source": [
    "## Chain-rule verification\n",
    "\n",
    "Finally we regard a Cartesian function of three\n",
    "variables as a function of the new coordinates.  Expanding the\n",
    "displayed coordinate formula should cancel every mixed derivative\n",
    "and leave the Cartesian Laplacian.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9eff2366",
   "metadata": {},
   "outputs": [],
   "source": [
    "def X : MathValue := r * sin θ * cos φ\n",
    "def Y : MathValue := r * sin θ * sin φ\n",
    "def Z : MathValue := r * cos θ\n",
    "def u : MathValue := function (X, Y, Z)\n",
    "\n",
    "def uR : MathValue := ∂/∂ u r\n",
    "def uRR : MathValue := ∂/∂ (∂/∂ u r) r\n",
    "def uΘ : MathValue := ∂/∂ u θ\n",
    "def uΘΘ : MathValue := ∂/∂ (∂/∂ u θ) θ\n",
    "def uΦΦ : MathValue := ∂/∂ (∂/∂ u φ) φ\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0d9433e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "$\\frac{\\partial^2 u}{\\partial 3^2} + \\frac{\\partial^2 u}{\\partial 2^2} + \\frac{\\partial^2 u}{\\partial 1^2}$"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "uRR + 2 * uR / r + uΘΘ / r^2\n",
    "  + cos θ * uΘ / (r^2 * sin θ)\n",
    "  + uΦΦ / (r^2 * (sin θ)^2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e95589e1",
   "metadata": {},
   "source": [
    "The expanded result is $u_{XX}+u_{YY}+u_{ZZ}$.  The first-derivative\n",
    "terms are exactly the corrections produced by the changing radial\n",
    "and angular scale factors.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Egison",
   "language": "egison",
   "name": "egison"
  },
  "language_info": {
   "codemirror_mode": "egison",
   "file_extension": ".egi",
   "mimetype": "text/x-egison",
   "name": "egison"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}