How to compute a gradient, a divergence or a curl?

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:

In [1]:
version()
Out[1]:
'SageMath version 9.0, Release Date: 2020-01-01'

First we set up the notebook to display math formulas using LaTeX formatting:

In [2]:
%display latex

First stage: introduce the Euclidean 3-space

Before evaluating some vector-field operators, one needs to define the arena in which vector fields live, namely the 3-dimensional Euclidean space $\mathbb{E}^3$. In SageMath, we declare it, along with the standard Cartesian coordinates $(x,y,z)$, as follows:

In [3]:
E.<x,y,z> = EuclideanSpace()
print(E)
E
Euclidean space E^3
Out[3]:

Thanks to the notation <x,y,z> in the above declaration, the coordinates $(x,y,z)$ are immediately available as three symbolic variables x, y and z (there is no need to type x, y, z = var('x y z')):

In [4]:
x is E.cartesian_coordinates()[1]
Out[4]:
In [5]:
y is E.cartesian_coordinates()[2]
Out[5]:
In [6]:
z is E.cartesian_coordinates()[3]
Out[6]:
In [7]:
type(z)
Out[7]:

Besides, $\mathbb{E}^3$ is endowed with the orthonormal vector frame $(e_x, e_y, e_z)$ associated with Cartesian coordinates:

In [8]:
E.frames()
Out[8]:

At this stage, this is the default vector frame on $\mathbb{E}^3$ (being the only vector frame introduced so far):

In [9]:
E.default_frame()
Out[9]:

Defining a vector field

We define a vector field on $\mathbb{E}^3$ from its components in the vector frame $(e_x,e_y,e_z)$:

In [10]:
v = E.vector_field(-y, x, sin(x*y*z), name='v')
v.display()
Out[10]:

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

In [11]:
v[1]
Out[11]:
In [12]:
v[:]
Out[12]:

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

In [13]:
p = E((3,-2,1), name='p')
print(p)
Point p on the Euclidean space E^3
In [14]:
p.coordinates()
Out[14]:
In [15]:
vp = v.at(p)
print(vp)
Vector v at Point p on the Euclidean space E^3
In [16]:
vp.display()
Out[16]:

Vector fields can be plotted (see the list of options for customizing the plot)

In [17]:
v.plot(max_range=1.5, scale=0.5)
Out[17]:

We may define a vector field with generic components $(u_x, u_y, u_z)$:

In [18]:
u = E.vector_field(function('u_x')(x,y,z),
                   function('u_y')(x,y,z),
                   function('u_z')(x,y,z),
                   name='u')
u.display()
Out[18]:
In [19]:
u[:]
Out[19]:

Its value at the point $p$ is then

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

How to compute various vector products?

Dot product

The dot (or scalar) product $u\cdot v$ of the vector fields $u$ and $v$ is obtained by the method dot_product, which admits dot as a shortcut alias:

In [21]:
u.dot(v) == u[1]*v[1] + u[2]*v[2] + u[3]*v[3]
Out[21]:

$s= u\cdot v$ is a scalar field, i.e. a map $\mathbb{E}^3 \to \mathbb{R}$:

In [22]:
s = u.dot(v)
s
Out[22]:
In [23]:
print(s)
Scalar field u.v on the Euclidean space E^3
In [24]:
s.display()
Out[24]:

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

In [25]:
s(p)
Out[25]:

Its coordinate expression is

In [26]:
s.expr()
Out[26]:

Norm

The norm $\|u\|$ of the vector field $u$ is defined in terms of the dot product by $\|u\| = \sqrt{u\cdot u}$:

In [27]:
norm(u) == sqrt(u.dot(u))
Out[27]:

It is a scalar field on $\mathbb{E}^3$:

In [28]:
s = norm(u)
print(s)
Scalar field |u| on the Euclidean space E^3
In [29]:
s.display()
Out[29]:
In [30]:
s.expr()
Out[30]:

For $v$, we have

In [31]:
norm(v).expr()
Out[31]:

Cross product

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

In [32]:
s = u.cross(v)
print(s)
Vector field u x v on the Euclidean space E^3
In [33]:
s.display()
Out[33]:

We can check the standard formulas expressing the cross product in terms of the components:

In [34]:
all([s[1] == u[2]*v[3] - u[3]*v[2],
     s[2] == u[3]*v[1] - u[1]*v[3],
     s[3] == u[1]*v[2] - u[2]*v[1]])
Out[34]:

Scalar triple product

Let us introduce a third vector field, $w$ say. 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 [35]:
w = E.vector_field(name='w')
w[1] = x*z
w[2] = y*z
w.display()
Out[35]:

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

In [36]:
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 [37]:
s.expr()
Out[37]:

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

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

How to evaluate the standard differential operators?

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

In [39]:
from sage.manifolds.operators import *

Gradient of a scalar field

We first introduce a scalar field, via its expression in terms of Cartesian coordinates; in this example, we consider some unspecified function of $(x,y,z)$:

In [40]:
F = E.scalar_field(function('f')(x,y,z), name='F')
F.display()
Out[40]:

The value of $F$ at a point:

In [41]:
F(p)
Out[41]:

The gradient of $F$:

In [42]:
print(grad(F))
Vector field grad(F) on the Euclidean space E^3
In [43]:
grad(F).display()
Out[43]:
In [44]:
norm(grad(F)).display()
Out[44]:

Divergence

The divergence of the vector field $u$:

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

For $v$ and $w$, we have

In [46]:
div(v).expr()
Out[46]:
In [47]:
div(w).expr()
Out[47]:

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

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

Curl

The curl of the vector field $u$:

In [49]:
s = curl(u)
print(s)
Vector field curl(u) on the Euclidean space E^3
In [50]:
s.display()
Out[50]:

To use the notation rot instead of curl, simply do

In [51]:
rot = curl

An alternative is

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

We have then

In [53]:
rot(u).display()
Out[53]:
In [54]:
rot(u) == curl(u)
Out[54]:

For $v$ and $w$, we have

In [55]:
curl(v).display()
Out[55]:
In [56]:
curl(w).display()
Out[56]:

The curl of a gradient is always zero:

In [57]:
curl(grad(F)).display()
Out[57]:

The divergence of a curl is always zero:

In [58]:
div(curl(u)).display()
Out[58]:

An identity valid for any scalar field $F$ and any vector field $u$ is $$ \mathrm{curl}(Fu) = \mathrm{grad}\, F\times u + F\, \mathrm{curl}\, u, $$ as we can check:

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

Laplacian

The Laplacian $\Delta F$ of a scalar field $F$ is another scalar field:

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

The following identity holds: $$\Delta F = \mathrm{div}\left(\mathrm{grad}\, F\right)$$ as we can check:

In [61]:
laplacian(F) == div(grad(F))
Out[61]:

The Laplacian $\Delta u$ of a vector field $u$ is another vector field:

In [62]:
Du = laplacian(u)
print(Du)
Vector field Delta(u) on the Euclidean space E^3

whose components are

In [63]:
Du.display()
Out[63]:

In the Cartesian frame, the components of $\Delta u$ are nothing but the (scalar) Laplacians of the components of $u$, as we can check:

In [64]:
e = E.cartesian_frame()
Du == sum(laplacian(u[[i]])*e[i] for i in E.irange())
Out[64]:

In the above formula, u[[i]] return the $i$-th component of $u$ as a scalar field, while u[i] would have returned the coordinate expression of this scalar field; besides, e is the Cartesian frame:

In [65]:
e[:]
Out[65]:

For the vector fields $v$ and $w$, we have

In [66]:
laplacian(v).display()
Out[66]:
In [67]:
laplacian(w).display()
Out[67]:

We have

In [68]:
curl(curl(u)).display()
Out[68]:
In [69]:
grad(div(u)).display()
Out[69]:

A famous identity is

$$ \mathrm{curl}\left(\mathrm{curl}\, u\right) = \mathrm{grad}\left(\mathrm{div}\, u\right) - \Delta u.$$

Let us check it:

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

How to customize various symbols?

Customizing the symbols of the orthonormal frame vectors

By default, the vectors of the orthonormal frame associated with Cartesian coordinates are denoted $(e_x,e_y,e_z)$:

In [71]:
frame = E.cartesian_frame()
frame
Out[71]:

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

In [72]:
frame.set_name('a', indices=('x', 'y', 'z'))
frame
Out[72]:
In [73]:
v.display()
Out[73]:
In [74]:
frame.set_name(('hx', 'hy', 'hz'), 
               latex_symbol=(r'\mathbf{\hat{x}}', r'\mathbf{\hat{y}}', 
                             r'\mathbf{\hat{z}}'))
frame
Out[74]:
In [75]:
v.display()
Out[75]:

Customizing the coordinate symbols

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

In [76]:
E.<x,y,z> = EuclideanSpace()

which resulted in the coordinate symbols $(x,y,z)$ and in the corresponding Python variables x, y and z (SageMath symbolic expressions). To use other symbols, for instance $(X,Y,Z)$, it suffices to create E as

In [77]:
E.<X,Y,Z> = EuclideanSpace()

We have then:

In [78]:
E.atlas()
Out[78]:
In [79]:
E.cartesian_frame()
Out[79]:
In [80]:
v = E.vector_field(-Y, X, sin(X*Y*Z), name='v')
v.display()
Out[80]:

By default the LaTeX symbols of the coordinate coincide with the letters given within the angle brackets. But this can be adjusted through the optional argument symbols of the function EuclideanSpace, which 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 [81]:
E.<xi,et,ze> = EuclideanSpace(symbols=r"xi:\xi et:\eta ze:\zeta")
E.atlas()
Out[81]:
In [82]:
v = E.vector_field(-et, xi, sin(xi*et*ze), name='v')
v.display()
Out[82]: