This notebook illustrates some vector calculus capabilities of SageMath within the 3-dimensional Euclidean space. The corresponding tools have been developed within the SageManifolds project.
Click here to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, via the command sage -n jupyter
NB: a version of SageMath at least equal to 8.3 is required to run this notebook:
version()
'SageMath version 8.8.beta2, Release Date: 2019-04-14'
First we set up the notebook to display math formulas using LaTeX formatting:
%display latex
We start by declaring the 3-dimensional Euclidean space $\mathbb{E}^3$, with $(r,\theta,\phi)$ as spherical coordinates:
E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
print(E)
E
Euclidean space E^3
$\mathbb{E}^3$ is endowed with the chart of spherical coordinates:
E.atlas()
as well as with the associated orthonormal vector frame $(e_r, e_\theta, e_\phi)$:
E.frames()
In the above output, $\left(\frac{\partial}{\partial r}, \frac{\partial}{\partial\theta}, \frac{\partial}{\partial \phi}\right)$ is the coordinate frame associated with $(r,\theta,\phi)$; it is not an orthonormal frame and will not be used below.
We define a vector field on $\mathbb{E}^3$ from its components in the orthonormal vector frame $(e_r,e_\theta,e_\phi)$:
v = E.vector_field(r*sin(2*ph)*sin(th)^2 + r,
r*sin(2*ph)*sin(th)*cos(th),
2*r*cos(ph)^2*sin(th), name='v')
v.display()
We can access to the components of $v$ via the square bracket operator:
v[1]
v[:]
A vector field can evaluated at any point of $\mathbb{E}^3$:
p = E((1, pi/2, pi), name='p')
print(p)
Point p on the Euclidean space E^3
p.coordinates()
vp = v.at(p)
print(vp)
Vector v at Point p on the Euclidean space E^3
vp.display()
We may define a vector field with generic components:
u = E.vector_field(function('u_r')(r,th,ph),
function('u_theta')(r,th,ph),
function('u_phi')(r,th,ph),
name='u')
u.display()
u[:]
up = u.at(p)
up.display()
s = u.dot(v)
s
print(s)
Scalar field u.v on the Euclidean space E^3
$s= u\cdot v$ is a scalar field, i.e. a map $\mathbb{E}^3 \rightarrow \mathbb{R}$:
s.display()
It maps points of $\mathbb{E}^3$ to real numbers:
s(p)
Its coordinate expression is
s.expr()
The norm of a vector field is
s = norm(u)
s
s.display()
s.expr()
The norm is related to the dot product by $\|u\|^2 = u\cdot u$, as we can check:
norm(u)^2 == u.dot(u)
For $v$, we have
norm(v).expr()
The cross product of $u$ by $v$ is obtained by the method cross_product
, which admits cross
as a shortcut alias:
s = u.cross(v)
print(s)
Vector field u x v on the Euclidean space E^3
s.display()
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:
w = E.vector_field(name='w')
w[1] = r
w.display()
The scalar triple product of the vector fields $u$, $v$ and $w$ is obtained as follows:
triple_product = E.scalar_triple_product()
s = triple_product(u, v, w)
print(s)
Scalar field epsilon(u,v,w) on the Euclidean space E^3
s.expr()
Let us check that the scalar triple product of $u$, $v$ and $w$ is $u\cdot(v\times w)$:
s == u.dot(v.cross(w))
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
(e.g. div(v)
):
from sage.manifolds.operators import *
We first introduce a scalar field, via its expression in terms of Cartesian coordinates; in this example, we consider a unspecified function of $(r,\theta,\phi)$:
F = E.scalar_field(function('f')(r,th,ph), name='F')
F.display()
The value of $F$ at a point:
F(p)
The gradient of $F$:
print(grad(F))
Vector field grad(F) on the Euclidean space E^3
grad(F).display()
norm(grad(F)).display()
The divergence of a vector field:
s = div(u)
s.display()
s.expr().expand()
For $v$ and $w$, we have
div(v).expr()
div(w).expr()
An identity valid for any scalar field $F$ and any vector field $u$:
div(F*u) == F*div(u) + u.dot(grad(F))
The curl of a vector field:
s = curl(u)
print(s)
Vector field curl(u) on the Euclidean space E^3
s.display()
To use the notation rot
instead of curl
, simply do
rot = curl
An alternative is
from sage.manifolds.operators import curl as rot
We have then
rot(u).display()
rot(u) == curl(u)
For $v$ and $w$, we have
curl(v).display()
curl(w).display()
The curl of a gradient is always zero:
curl(grad(F)).display()
The divergence of a curl is always zero:
div(curl(u)).display()
An identity valid for any scalar field $F$ and any vector field $u$:
curl(F*u) == grad(F).cross(u) + F*curl(u)
The Laplacian of a scalar field:
s = laplacian(F)
s.display()
s.expr().expand()
For a scalar field, the Laplacian is nothing but the divergence of the gradient:
laplacian(F) == div(grad(F))
The Laplacian of a vector field:
Du = laplacian(u)
Du.display()
Since this expression is quite lengthy, we may ask for a display component by component:
Du.display_comp()
We may expand each component:
for i in E.irange():
Du[i].expand()
Du.display_comp()
Du[1]
Du[2]
Du[3]
As a test, we may check that these formulas coincide with those of Wikipedia's article Del in cylindrical and spherical coordinates.
For $v$ and $w$, we have
laplacian(v).display()
laplacian(w).display()
We have
curl(curl(u)).display()
grad(div(u)).display()
and we may check a famous identity:
curl(curl(u)) == grad(div(u)) - laplacian(u)
frame = E.spherical_frame()
frame
But this can be changed, thanks to the method set_name
:
frame.set_name('a', indices=('r', 'th', 'ph'),
latex_indices=('r', r'\theta', r'\phi'))
frame
v.display()
frame.set_name(('hr', 'hth', 'hph'),
latex_symbol=(r'\hat{r}', r'\hat{\theta}', r'\hat{\phi}'))
frame
v.display()
The coordinates symbols are defined within the angle brackets <...>
at the construction of the Euclidean space. Above we did
E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
which resulted in the coordinate symbols $(r,\theta,\phi)$ and in the corresponding Python variables r
, th
and ph
(SageMath symbolic expressions). Using other symbols, for instance $(R,\Theta,\Phi)$, 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 (:
):
E.<R,Th,Ph> = EuclideanSpace(coordinates='spherical', symbols=r'R Th:\Theta Ph:\Phi')
We have then
E.atlas()
E.frames()
E.spherical_frame()
v = E.vector_field(R*sin(2*Ph)*sin(Th)^2 + R,
R*sin(2*Ph)*sin(Th)*cos(Th),
2*R*cos(Ph)^2*sin(Th), name='v')
v.display()