Vector calculus with SageMath¶

2. Using spherical coordinates¶

This notebook illustrates operations on vector fields on Euclidean spaces, as introduced in Trac ticket #24623.

In [1]:
version()

Out[1]:
'SageMath version 8.2.rc4, Release Date: 2018-04-20'
In [2]:
%display latex

In [3]:
from sage.manifolds.operators import *   # to get the operators grad, div, curl, etc.


The 3-dimensional Euclidean space¶

We start by declaring the 3-dimensional Euclidean space $\mathbb{E}^3$, with $(r,\theta,\phi)$ as spherical coordinates:

In [4]:
E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
print(E)
E

Euclidean space E^3

Out[4]:

$\mathbb{E}^3$ is endowed with the chart of spherical coordinates:

In [5]:
E.atlas()

Out[5]:

as well as with the associated orthonormal vector frame $(e_r, e_\theta, e_\phi)$:

In [6]:
E.frames()

Out[6]:

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.

Vector fields¶

We define a vector field on $\mathbb{E}^3$ from its components in the orthonormal vector frame $(e_r,e_\theta,e_\phi)$:

In [7]:
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()

Out[7]:

We can access to the components of $v$ via the square bracket operator:

In [8]:
v[1]

Out[8]:
In [9]:
v[:]

Out[9]:

A vector field can evaluated at any point of $\mathbb{E}^3$:

In [10]:
p = E((1, pi/2, pi), name='p')
print(p)

Point p on the Euclidean space E^3

In [11]:
p.coordinates()

Out[11]:
In [12]:
vp = v.at(p)
print(vp)

Vector v at Point p on the Euclidean space E^3

In [13]:
vp.display()

Out[13]:

We may define a vector field with generic components:

In [14]:
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()

Out[14]:
In [15]:
u[:]

Out[15]:
In [16]:
up = u.at(p)
up.display()

Out[16]:

Algebraic operations on vector fields¶

Dot product¶

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:

In [17]:
s = u.dot(v)
s

Out[17]:
In [18]:
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}$:

In [19]:
s.display()

Out[19]:

It maps point of $\mathbb{E}^3$ to real numbers:

In [20]:
s(p)

Out[20]:

Its coordinate expression is

In [21]:
s.expr()

Out[21]:

Norm¶

The norm of a vector field is

In [22]:
s = norm(u)
s

Out[22]:
In [23]:
s.display()

Out[23]:
In [24]:
s.expr()

Out[24]:

The norm is related to the dot product by $\|u\|^2 = u\cdot u$, as we can check:

In [25]:
norm(u)^2 == u.dot(u)

Out[25]:

For $v$, we have:

In [26]:
norm(v).expr()

Out[26]:

Cross product¶

The cross product of $u$ by $v$ is obtained by the method cross_product, which admits cross as a shortcut alias:

In [27]:
s = u.cross(v)
print(s)

Vector field u x v on the Euclidean space E^3

In [28]:
s.display()

Out[28]:

Scalar triple product¶

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:

In [29]:
w = E.vector_field(name='w')
w[1] = r
w.display()

Out[29]:

The scalar triple product of the vector fields $u$, $v$ and $w$ is obtained as follows:

In [30]:
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

In [31]:
s.expr()

Out[31]:

Let us check that the scalar triple product of $u$, $v$ and $w$ is $u\cdot(v\times w)$:

In [32]:
s == u.dot(v.cross(w))

Out[32]:

Differential operators¶

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)$:

In [33]:
F = E.scalar_field(function('f')(r,th,ph), name='F')
F.display()

Out[33]:

The value of $F$ at a point:

In [34]:
F(p)

Out[34]:

The gradient of $F$:

In [35]:
print(grad(F))

Vector field grad(F) on the Euclidean space E^3

In [36]:
grad(F).display()

Out[36]:
In [37]:
norm(grad(F)).display()

Out[37]:

Divergence¶

The divergence of a vector field:

In [38]:
s = div(u)
s.display()

Out[38]:
In [39]:
s.expr().expand()

Out[39]:

For $v$ and $w$, we have

In [40]:
div(v).expr()

Out[40]:
In [41]:
div(w).expr()

Out[41]:

An identity valid for any scalar field $F$ and any vector field $u$:

In [42]:
div(F*u) == F*div(u) + u.dot(grad(F))

Out[42]:

Curl¶

The curl of a vector field:

In [43]:
s = curl(u)
print(s)

Vector field curl(u) on the Euclidean space E^3

In [44]:
s.display()

Out[44]:

To use the notation rot instead of curl, simply do

In [45]:
rot = curl


An alternative is

In [46]:
from sage.manifolds.operators import curl as rot


We have then

In [47]:
rot(u).display()

Out[47]:
In [48]:
rot(u) == curl(u)

Out[48]:

For $v$ and $w$, we have:

In [49]:
curl(v).display()

Out[49]:
In [50]:
curl(w).display()

Out[50]:

The curl of a gradient is always zero:

In [51]:
curl(grad(F)).display()

Out[51]:

The divergence of a curl is always zero:

In [52]:
div(curl(u)).display()

Out[52]:

An identity valid for any scalar field $F$ and any vector field $u$:

In [53]:
curl(F*u) == grad(F).cross(u) + F*curl(u)

Out[53]:

Laplacian¶

The Laplacian of a scalar field:

In [54]:
s = laplacian(F)
s.display()

Out[54]:
In [55]:
s.expr().expand()

Out[55]:

For a scalar field, the Laplacian is nothing but the divergence of the gradient:

In [56]:
laplacian(F) == div(grad(F))

Out[56]:

The Laplacian of a vector field:

In [57]:
Du = laplacian(u)
Du.display()

Out[57]:

Since this expression is quite lengthy, we may ask for a display component by component:

In [58]:
Du.display_comp()

Out[58]:

We may expand each component:

In [59]:
for i in E.irange():
Du[i].expand()
Du.display_comp()

Out[59]:
In [60]:
Du[1]

Out[60]:
In [61]:
Du[2]

Out[61]:
In [62]:
Du[3]

Out[62]:

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

In [63]:
laplacian(v).display()

Out[63]:
In [64]:
laplacian(w).display()

Out[64]:

We have:

In [65]:
curl(curl(u)).display()

Out[65]:
In [66]:
grad(div(u)).display()

Out[66]:

and we may check a famous identity:

In [67]:
curl(curl(u)) == grad(div(u)) - laplacian(u)

Out[67]:

Customizations¶

Customizing the symbols of the orthonormal frame vectors¶

By default, the vectors of the orthonormal frame associated with spherical coordinates are denoted $(e_r,e_\theta,e_\phi)$:

In [68]:
frame = E.spherical_frame()
frame

Out[68]:

But this can be changed, thanks to the method set_name:

In [69]:
frame.set_name('a', indices=('r', 'th', 'ph'),
latex_indices=('r', r'\theta', r'\phi'))
frame

Out[69]:
In [70]:
v.display()

Out[70]:
In [71]:
frame.set_name(('hr', 'hth', 'hph'),
latex_symbol=(r'\hat{r}', r'\hat{\theta}', r'\hat{\phi}'))
frame

Out[71]:
In [72]:
v.display()

Out[72]:

Customizing the coordinate symbols¶

The coordinates symbols are defined within the angle brackets <...> at the construction of the Euclidean space. Above we did

In [73]:
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 (:):

In [74]:
E.<R,Th,Ph> = EuclideanSpace(coordinates='spherical', symbols=r'R Th:\Theta Ph:\Phi')


We have then

In [75]:
E.atlas()

Out[75]:
In [76]:
E.frames()

Out[76]:
In [77]:
E.spherical_frame()

Out[77]:
In [78]:
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()

Out[78]:
In [ ]: