{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Vector calculus with SageMath\n",
"\n",
"## Part 3: Using cylindrical coordinates\n",
"\n",
"This notebook illustrates some vector calculus capabilities of SageMath within the 3-dimensional Euclidean space. The corresponding tools have been developed within\n",
"the [SageManifolds](https://sagemanifolds.obspm.fr) project.\n",
"\n",
"Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Notebooks/SM_vector_calc_cylindrical.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, via the command `sage -n jupyter`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*NB:* a version of SageMath at least equal to 8.3 is required to run this notebook:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'SageMath version 8.8.beta2, Release Date: 2019-04-14'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"version()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we set up the notebook to display math formulas using LaTeX formatting:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%display latex"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The 3-dimensional Euclidean space\n",
"\n",
"We start by declaring the 3-dimensional Euclidean space $\\mathbb{E}^3$, with $(\\rho,\\phi,z)$ as cylindrical coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Euclidean space E^3\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Euclidean space E^3"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace(coordinates='cylindrical')\n",
"print(E)\n",
"E"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\mathbb{E}^3$ is endowed with the chart of cylindrical coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[Chart (E^3, (rh, ph, z))]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E.atlas()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"as well as with the associated orthonormal vector frame $(e_\\rho, e_\\phi, e_z)$:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[Coordinate frame (E^3, (d/drh,d/dph,d/dz)),\n",
" Vector frame (E^3, (e_rh,e_ph,e_z))]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E.frames()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the above output, $\\left(\\frac{\\partial}{\\partial\\rho}, \\frac{\\partial}{\\partial\\phi}, \\frac{\\partial}{\\partial z}\\right)$ is the coordinate frame associated with $(\\rho,\\phi,z)$; it is not an orthonormal frame and will not be used below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Vector fields"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We define a vector field on $\\mathbb{E}^3$ from its components in the orthonormal vector frame $(e_\\rho,e_\\phi,e_z)$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"v = rh*(sin(2*ph) + 1) e_rh + 2*rh*cos(ph)^2 e_ph + z e_z"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v = E.vector_field(rh*(1+sin(2*ph)), 2*rh*cos(ph)^2, z,\n",
" name='v')\n",
"v.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can access to the components of $v$ via the square bracket operator:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"rh*(sin(2*ph) + 1)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v[1]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[rh*(sin(2*ph) + 1), 2*rh*cos(ph)^2, z]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v[:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A vector field can evaluated at any point of $\\mathbb{E}^3$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Point p on the Euclidean space E^3\n"
]
}
],
"source": [
"p = E((1, pi, 0), name='p')\n",
"print(p)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(1, pi, 0)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p.coordinates()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Vector v at Point p on the Euclidean space E^3\n"
]
}
],
"source": [
"vp = v.at(p)\n",
"print(vp)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"v = e_rh + 2 e_ph"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vp.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We may define a vector field with generic components:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"u = u_rho(rh, ph, z) e_rh + u_phi(rh, ph, z) e_ph + u_z(rh, ph, z) e_z"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = E.vector_field(function('u_rho')(rh,ph,z),\n",
" function('u_phi')(rh,ph,z),\n",
" function('u_z')(rh,ph,z),\n",
" name='u')\n",
"u.display()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[u_rho(rh, ph, z), u_phi(rh, ph, z), u_z(rh, ph, z)]"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u[:]"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"u = u_rho(1, pi, 0) e_rh + u_phi(1, pi, 0) e_ph + u_z(1, pi, 0) e_z"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"up = u.at(p)\n",
"up.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Algebraic operations on vector fields\n",
"\n",
"### Dot product\n",
"\n",
"The dot (or scalar) product of the vector fields $u$ and $v$ is obtained by the method `dot_product`, which admits `dot` as a shortcut alias:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Scalar field u.v on the Euclidean space E^3"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = u.dot(v)\n",
"s"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Scalar field u.v on the Euclidean space E^3\n"
]
}
],
"source": [
"print(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$s= u\\cdot v$ is a *scalar field*, i.e. a map $\\mathbb{E}^3 \\rightarrow \\mathbb{R}$:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"u.v: E^3 --> R\n",
" (rh, ph, z) |--> 2*rh*cos(ph)^2*u_phi(rh, ph, z) + (2*cos(ph)*sin(ph) + 1)*rh*u_rho(rh, ph, z) + z*u_z(rh, ph, z)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It maps points of $\\mathbb{E}^3$ to real numbers:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"2*u_phi(1, pi, 0) + u_rho(1, pi, 0)"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Its coordinate expression is"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"2*rh*cos(ph)^2*u_phi(rh, ph, z) + (2*cos(ph)*sin(ph) + 1)*rh*u_rho(rh, ph, z) + z*u_z(rh, ph, z)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Norm\n",
"\n",
"The norm of a vector field is "
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Scalar field |u| on the Euclidean space E^3"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = norm(u)\n",
"s"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"|u|: E^3 --> R\n",
" (rh, ph, z) |--> sqrt(u_phi(rh, ph, z)^2 + u_rho(rh, ph, z)^2 + u_z(rh, ph, z)^2)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.display()"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"sqrt(u_phi(rh, ph, z)^2 + u_rho(rh, ph, z)^2 + u_z(rh, ph, z)^2)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The norm is related to the dot product by $\\|u\\|^2 = u\\cdot u$, as we can check:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(u)^2 == u.dot(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $v$, we have"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"sqrt((4*cos(ph)^2 + 4*cos(ph)*sin(ph) + 1)*rh^2 + z^2)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(v).expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cross product\n",
"\n",
"The cross product of $u$ by $v$ is obtained by the method `cross_product`, which admits `cross` as a shortcut alias:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Vector field u x v on the Euclidean space E^3\n"
]
}
],
"source": [
"s = u.cross(v)\n",
"print(s)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"u x v = (-2*rh*cos(ph)^2*u_z(rh, ph, z) + z*u_phi(rh, ph, z)) e_rh + ((2*cos(ph)*sin(ph) + 1)*rh*u_z(rh, ph, z) - z*u_rho(rh, ph, z)) e_ph + (2*rh*cos(ph)^2*u_rho(rh, ph, z) - (2*cos(ph)*sin(ph) + 1)*rh*u_phi(rh, ph, z)) e_z"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Scalar triple product\n",
"\n",
"Let us introduce a third vector field. As a example, we do not pass the components as arguments of `vector_field`, as we did for $u$ and $v$; instead, we set them in a second stage, via the square bracket operator, any unset component being assumed to be zero:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"w = rh e_rh + z e_z"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w = E.vector_field(name='w')\n",
"w[1] = rh\n",
"w[3] = z\n",
"w.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The scalar triple product of the vector fields $u$, $v$ and $w$ is obtained as follows: "
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Scalar field epsilon(u,v,w) on the Euclidean space E^3\n"
]
}
],
"source": [
"triple_product = E.scalar_triple_product()\n",
"s = triple_product(u, v, w)\n",
"print(s)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"-2*rh^2*cos(ph)^2*u_z(rh, ph, z) - 2*(rh*cos(ph)*sin(ph)*u_phi(rh, ph, z) - rh*cos(ph)^2*u_rho(rh, ph, z))*z"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us check that the scalar triple product of $u$, $v$ and $w$ is $u\\cdot(v\\times w)$:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s == u.dot(v.cross(w))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Differential operators"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While the standard operators $\\mathrm{grad}$, $\\mathrm{div}$, $\\mathrm{curl}$, etc. involved in vector calculus are accessible via the dot notation (e.g. `v.div()`), let us import functions `grad`, `div`, `curl`, etc. that allow for using standard mathematical notations\n",
"(e.g. `div(v)`):"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"from sage.manifolds.operators import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gradient of a scalar field\n",
"\n",
"We first introduce a scalar field, via its expression in terms of Cartesian coordinates; in this example, we consider a unspecified function of $(\\rho,\\phi,z)$:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"F: E^3 --> R\n",
" (rh, ph, z) |--> f(rh, ph, z)"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"F = E.scalar_field(function('f')(rh,ph,z), name='F')\n",
"F.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The value of $F$ at a point:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"f(1, pi, 0)"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"F(p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The gradient of $F$:"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Vector field grad(F) on the Euclidean space E^3\n"
]
}
],
"source": [
"print(grad(F))"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"grad(F) = d(f)/drh e_rh + d(f)/dph/rh e_ph + d(f)/dz e_z"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"grad(F).display()"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"|grad(F)|: E^3 --> R\n",
" (rh, ph, z) |--> sqrt(rh^2*(d(f)/drh)^2 + rh^2*(d(f)/dz)^2 + (d(f)/dph)^2)/rh"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(grad(F)).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Divergence \n",
"\n",
"The divergence of a vector field:"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"div(u): E^3 --> R\n",
" (rh, ph, z) |--> (rh*d(u_rho)/drh + rh*d(u_z)/dz + u_rho(rh, ph, z) + d(u_phi)/dph)/rh"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = div(u)\n",
"s.display()"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"u_rho(rh, ph, z)/rh + diff(u_phi(rh, ph, z), ph)/rh + diff(u_rho(rh, ph, z), rh) + diff(u_z(rh, ph, z), z)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.expr().expand()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $v$ and $w$, we have"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"3"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"div(v).expr()"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"3"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"div(w).expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An identity valid for any scalar field $F$ and any vector field $u$:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"div(F*u) == F*div(u) + u.dot(grad(F))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Curl\n",
"\n",
"The curl of a vector field:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Vector field curl(u) on the Euclidean space E^3\n"
]
}
],
"source": [
"s = curl(u)\n",
"print(s)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(u) = -(rh*d(u_phi)/dz - d(u_z)/dph)/rh e_rh + (d(u_rho)/dz - d(u_z)/drh) e_ph + (rh*d(u_phi)/drh + u_phi(rh, ph, z) - d(u_rho)/dph)/rh e_z"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To use the notation `rot` instead of `curl`, simply do"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"rot = curl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative is"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"from sage.manifolds.operators import curl as rot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have then"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(u) = -(rh*d(u_phi)/dz - d(u_z)/dph)/rh e_rh + (d(u_rho)/dz - d(u_z)/drh) e_ph + (rh*d(u_phi)/drh + u_phi(rh, ph, z) - d(u_rho)/dph)/rh e_z"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rot(u).display()"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rot(u) == curl(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $v$ and $w$, we have"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(v) = 2 e_z"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(v).display()"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(w) = 0"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(w).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The curl of a gradient is always zero:"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(grad(F)) = 0"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(grad(F)).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The divergence of a curl is always zero:"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"div(curl(u)): E^3 --> R\n",
" (rh, ph, z) |--> 0"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"div(curl(u)).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An identity valid for any scalar field $F$ and any vector field $u$:"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(F*u) == grad(F).cross(u) + F*curl(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Laplacian"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Laplacian of a scalar field:"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(F): E^3 --> R\n",
" (rh, ph, z) |--> (rh^2*d^2(f)/drh^2 + rh^2*d^2(f)/dz^2 + rh*d(f)/drh + d^2(f)/dph^2)/rh^2"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = laplacian(F)\n",
"s.display()"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"diff(f(rh, ph, z), rh)/rh + diff(f(rh, ph, z), ph, ph)/rh^2 + diff(f(rh, ph, z), rh, rh) + diff(f(rh, ph, z), z, z)"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.expr().expand()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a scalar field, the Laplacian is nothing but the divergence of the gradient:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"laplacian(F) == div(grad(F))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Laplacian of a vector field:"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(u) = (rh^2*d^2(u_rho)/drh^2 + rh^2*d^2(u_rho)/dz^2 + rh*d(u_rho)/drh - u_rho(rh, ph, z) - 2*d(u_phi)/dph + d^2(u_rho)/dph^2)/rh^2 e_rh + (rh^2*d^2(u_phi)/drh^2 + rh^2*d^2(u_phi)/dz^2 + rh*d(u_phi)/drh - u_phi(rh, ph, z) + d^2(u_phi)/dph^2 + 2*d(u_rho)/dph)/rh^2 e_ph + (rh^2*d^2(u_z)/drh^2 + rh^2*d^2(u_z)/dz^2 + rh*d(u_z)/drh + d^2(u_z)/dph^2)/rh^2 e_z"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Du = laplacian(u)\n",
"Du.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since this expression is quite lengthy, we may ask for a display component by component:"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(u)^1 = (rh^2*d^2(u_rho)/drh^2 + rh^2*d^2(u_rho)/dz^2 + rh*d(u_rho)/drh - u_rho(rh, ph, z) - 2*d(u_phi)/dph + d^2(u_rho)/dph^2)/rh^2 \n",
"Delta(u)^2 = (rh^2*d^2(u_phi)/drh^2 + rh^2*d^2(u_phi)/dz^2 + rh*d(u_phi)/drh - u_phi(rh, ph, z) + d^2(u_phi)/dph^2 + 2*d(u_rho)/dph)/rh^2 \n",
"Delta(u)^3 = (rh^2*d^2(u_z)/drh^2 + rh^2*d^2(u_z)/dz^2 + rh*d(u_z)/drh + d^2(u_z)/dph^2)/rh^2 "
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Du.display_comp()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We may expand each component:"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(u)^1 = d(u_rho)/drh/rh - u_rho(rh, ph, z)/rh^2 - 2*d(u_phi)/dph/rh^2 + d^2(u_rho)/dph^2/rh^2 + d^2(u_rho)/drh^2 + d^2(u_rho)/dz^2 \n",
"Delta(u)^2 = d(u_phi)/drh/rh - u_phi(rh, ph, z)/rh^2 + d^2(u_phi)/dph^2/rh^2 + 2*d(u_rho)/dph/rh^2 + d^2(u_phi)/drh^2 + d^2(u_phi)/dz^2 \n",
"Delta(u)^3 = d(u_z)/drh/rh + d^2(u_z)/dph^2/rh^2 + d^2(u_z)/drh^2 + d^2(u_z)/dz^2 "
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"for i in E.irange():\n",
" Du[i].expand()\n",
"Du.display_comp()"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"d(u_rho)/drh/rh - u_rho(rh, ph, z)/rh^2 - 2*d(u_phi)/dph/rh^2 + d^2(u_rho)/dph^2/rh^2 + d^2(u_rho)/drh^2 + d^2(u_rho)/dz^2"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Du[1]"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"d(u_phi)/drh/rh - u_phi(rh, ph, z)/rh^2 + d^2(u_phi)/dph^2/rh^2 + 2*d(u_rho)/dph/rh^2 + d^2(u_phi)/drh^2 + d^2(u_phi)/dz^2"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Du[2]"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"d(u_z)/drh/rh + d^2(u_z)/dph^2/rh^2 + d^2(u_z)/drh^2 + d^2(u_z)/dz^2"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Du[3]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a test, we may check that these formulas coincide with those of\n",
"Wikipedia's article [*Del in cylindrical and spherical coordinates*](https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates#Del_formula)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $v$ and $w$, we have "
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(v) = 0"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"laplacian(v).display()"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Delta(w) = 0"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"laplacian(w).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"curl(curl(u)) = -(rh^2*d^2(u_rho)/dz^2 - rh^2*d^2(u_z)/drhdz - rh*d^2(u_phi)/drhdph - d(u_phi)/dph + d^2(u_rho)/dph^2)/rh^2 e_rh - (rh^2*d^2(u_phi)/drh^2 + rh^2*d^2(u_phi)/dz^2 + rh*d(u_phi)/drh - rh*d^2(u_rho)/drhdph - rh*d^2(u_z)/dphdz - u_phi(rh, ph, z) + d(u_rho)/dph)/rh^2 e_ph + (rh^2*d^2(u_rho)/drhdz - rh^2*d^2(u_z)/drh^2 + rh*d^2(u_phi)/dphdz + rh*d(u_rho)/dz - rh*d(u_z)/drh - d^2(u_z)/dph^2)/rh^2 e_z"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(curl(u)).display()"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"grad(div(u)) = (rh^2*d^2(u_rho)/drh^2 + rh^2*d^2(u_z)/drhdz + rh*d^2(u_phi)/drhdph + rh*d(u_rho)/drh - u_rho(rh, ph, z) - d(u_phi)/dph)/rh^2 e_rh + (rh*d^2(u_rho)/drhdph + rh*d^2(u_z)/dphdz + d^2(u_phi)/dph^2 + d(u_rho)/dph)/rh^2 e_ph + (rh*d^2(u_rho)/drhdz + rh*d^2(u_z)/dz^2 + d^2(u_phi)/dphdz + d(u_rho)/dz)/rh e_z"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"grad(div(u)).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and we may check a famous identity:"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"curl(curl(u)) == grad(div(u)) - laplacian(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Customizations\n",
"\n",
"### Customizing the symbols of the orthonormal frame vectors\n",
"\n",
"By default, the vectors of the orthonormal frame associated with cylindrical coordinates are denoted $(e_\\rho,e_\\phi,z)$:"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Vector frame (E^3, (e_rh,e_ph,e_z))"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"frame = E.cylindrical_frame()\n",
"frame"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But this can be changed, thanks to the method `set_name`:"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Vector frame (E^3, (a_rh,a_ph,a_z))"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"frame.set_name('a', indices=('rh', 'ph', 'z'),\n",
" latex_indices=(r'\\rho', r'\\phi', 'z'))\n",
"frame"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"v = rh*(sin(2*ph) + 1) a_rh + 2*rh*cos(ph)^2 a_ph + z a_z"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v.display()"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Vector frame (E^3, (hrh,hph,hz))"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"frame.set_name(('hrh', 'hph', 'hz'), \n",
" latex_symbol=(r'\\hat{\\rho}', r'\\hat{\\phi}', r'\\hat{z}'))\n",
"frame"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"v = rh*(sin(2*ph) + 1) hrh + 2*rh*cos(ph)^2 hph + z hz"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Customizing the coordinate symbols\n",
"\n",
"The coordinates symbols are defined within the angle brackets `<...>` at the construction of the Euclidean space. Above we did "
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [],
"source": [
"E. = EuclideanSpace(coordinates='cylindrical')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which resulted in the coordinate symbols $(\\rho,\\phi,z)$ and in the corresponding Python variables `rh`, `ph` and `z` (SageMath symbolic expressions). Using other symbols, for instance $(R,\\Phi,Z)$, is possible through the optional argument `symbols` of the function `EuclideanSpace`. It has to be a string, usually prefixed by `r` (for *raw* string, in order to allow for the backslash character of LaTeX expressions). This string contains the coordinate fields separated by a blank space; each field contains the coordinate’s text symbol and possibly the coordinate’s LaTeX symbol (when the latter is different from the text symbol), both symbols being separated by a colon (`:`):"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [],
"source": [
"E. = EuclideanSpace(coordinates='cylindrical', symbols=r'R Ph:\\Phi Z')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have then"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[Chart (E^3, (R, Ph, Z))]"
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E.atlas()"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[Coordinate frame (E^3, (d/dR,d/dPh,d/dZ)), Vector frame (E^3, (e_R,e_Ph,e_Z))]"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E.frames()"
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Vector frame (E^3, (e_R,e_Ph,e_Z))"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E.cylindrical_frame()"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"v = R*(sin(2*Ph) + 1) e_R + 2*R*cos(Ph)^2 e_Ph + Z e_Z"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v = E.vector_field(R*(1+sin(2*Ph)), 2*R*cos(Ph)^2, Z,\n",
" name='v')\n",
"v.display()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 8.8.beta2",
"language": "",
"name": "sagemath"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}