version()
'SageMath version 8.2.rc4, Release Date: 2018-04-20'
%display latex
from sage.manifolds.operators import * # to get the operators grad, div, curl, etc.
We start by declaring the 3-dimensional Euclidean space $\mathbb{E}^3$, with $(\rho,\phi,z)$ as spherical coordinates:
E.<rh,ph,z> = EuclideanSpace(coordinates='cylindrical')
print(E)
E
Euclidean space E^3
$\mathbb{E}^3$ is endowed with the chart of cylindrical coordinates:
E.atlas()
as well as with the associated orthonormal vector frame $(e_\rho, e_\phi, e_z)$:
E.frames()
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.
We define a vector field on $\mathbb{E}^3$ from its components in the orthonormal vector frame $(e_\rho,e_\phi,e_z)$:
v = E.vector_field(rh*(1+sin(2*ph)), 2*rh*cos(ph)^2, z,
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, 0), 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_rho')(rh,ph,z),
function('u_phi')(rh,ph,z),
function('u_z')(rh,ph,z),
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 point 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] = rh
w[3] = z
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))
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)$:
F = E.scalar_field(function('f')(rh,ph,z), 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.cylindrical_frame()
frame
But this can be changed, thanks to the method set_name
:
frame.set_name('a', indices=('rh', 'ph', 'z'),
latex_indices=(r'\rho', r'\phi', 'z'))
frame
v.display()
frame.set_name(('hrh', 'hph', 'hz'),
latex_symbol=(r'\hat{\rho}', r'\hat{\phi}', r'\hat{z}'))
frame
v.display()
The coordinates symbols are defined within the angle brackets <...>
at the construction of the Euclidean space. Above we did
E.<rh,ph,z> = EuclideanSpace(coordinates='cylindrical')
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 (:
):
E.<R,Ph,Z> = EuclideanSpace(coordinates='cylindrical', symbols=r'R Ph:\Phi Z')
We have then
E.atlas()
E.frames()
E.cylindrical_frame()
v = E.vector_field(R*(1+sin(2*Ph)), 2*R*cos(Ph)^2, Z,
name='v')
v.display()