How to compute a gradient, a divergence or a curl¶
This tutorial introduces some vector calculus capabilities of SageMath within the 3-dimensional Euclidean space. The corresponding tools have been developed via the SageManifolds project.
The tutorial is also available as a Jupyter notebook, either
passive (nbviewer)
or interactive (binder).
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)\), via EuclideanSpace:
sage: E.<x,y,z> = EuclideanSpace()
sage: E
Euclidean space E^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 declare them via var(), i.e. to type
x, y, z = var('x y z')):
sage: x is E.cartesian_coordinates()[1]
True
sage: y is E.cartesian_coordinates()[2]
True
sage: z is E.cartesian_coordinates()[3]
True
sage: type(z)
<type 'sage.symbolic.expression.Expression'>
Besides, \(\mathbb{E}^3\) is endowed with the orthonormal vector frame \((e_x, e_y, e_z)\) associated with Cartesian coordinates:
sage: E.frames()
[Coordinate frame (E^3, (e_x,e_y,e_z))]
At this stage, this is the default vector frame on \(\mathbb{E}^3\) (being the only vector frame introduced so far):
sage: E.default_frame()
Coordinate frame (E^3, (e_x,e_y,e_z))
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)\):
sage: v = E.vector_field(-y, x, sin(x*y*z), name='v')
sage: v.display()
v = -y e_x + x e_y + sin(x*y*z) e_z
We can access to the components of \(v\) via the square bracket operator:
sage: v[1]
-y
sage: v[:]
[-y, x, sin(x*y*z)]
A vector field can evaluated at any point of \(\mathbb{E}^3\):
sage: p = E((3,-2,1), name='p')
sage: p
Point p on the Euclidean space E^3
sage: p.coordinates()
(3, -2, 1)
sage: vp = v.at(p)
sage: vp
Vector v at Point p on the Euclidean space E^3
sage: vp.display()
v = 2 e_x + 3 e_y - sin(6) e_z
Vector fields can be plotted:
sage: v.plot(max_range=1.5, scale=0.5)
Graphics3d Object
For customizing the plot, see the list of options in the documentation of
plot().
For instance, to get a view of the orthogonal projection of \(v\) in the plane
\(y=1\), do:
sage: v.plot(fixed_coords={y: 1}, ambient_coords=(x,z), max_range=1.5,
....: scale=0.25)
Graphics object consisting of 81 graphics primitives
We may define a vector field \(u\) with generic components \((u_x, u_y, y_z)\):
sage: u = E.vector_field(function('u_x')(x,y,z),
....: function('u_y')(x,y,z),
....: function('u_z')(x,y,z),
....: name='u')
sage: u.display()
u = u_x(x, y, z) e_x + u_y(x, y, z) e_y + u_z(x, y, z) e_z
sage: u[:]
[u_x(x, y, z), u_y(x, y, z), u_z(x, y, z)]
Its value at the point \(p\) is then:
sage: up = u.at(p)
sage: up.display()
u = u_x(3, -2, 1) e_x + u_y(3, -2, 1) e_y + u_z(3, -2, 1) e_z
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:
sage: u.dot(v) == u[1]*v[1] + u[2]*v[2] + u[3]*v[3]
True
\(s= u\cdot v\) is a scalar field, i.e. a map \(\mathbb{E}^3 \to \mathbb{R}\):
sage: s = u.dot(v)
sage: s
Scalar field u.v on the Euclidean space E^3
sage: s.display()
u.v: E^3 → ℝ
(x, y, z) ↦ -y*u_x(x, y, z) + x*u_y(x, y, z) + sin(x*y*z)*u_z(x, y, z)
It maps points of \(\mathbb{E}^3\) to real numbers:
sage: s(p)
-sin(6)*u_z(3, -2, 1) + 2*u_x(3, -2, 1) + 3*u_y(3, -2, 1)
Its coordinate expression is:
sage: s.expr()
-y*u_x(x, y, z) + x*u_y(x, y, z) + sin(x*y*z)*u_z(x, y, z)
Norm¶
The norm \(\|u\|\) of the vector field \(u\) is defined in terms of the dot product by \(\|u\| = \sqrt{u\cdot u}\):
sage: norm(u) == sqrt(u.dot(u))
True
It is a scalar field on \(\mathbb{E}^3\):
sage: s = norm(u)
sage: s
Scalar field |u| on the Euclidean space E^3
sage: s.display()
|u|: E^3 → ℝ
(x, y, z) ↦ sqrt(u_x(x, y, z)^2 + u_y(x, y, z)^2 + u_z(x, y, z)^2)
sage: s.expr()
sqrt(u_x(x, y, z)^2 + u_y(x, y, z)^2 + u_z(x, y, z)^2)
For \(v\), we have:
sage: norm(v).expr()
sqrt(x^2 + y^2 + sin(x*y*z)^2)
Cross product¶
The cross product \(u\times v\) is obtained by the method
cross_product(),
which admits cross() as a shortcut alias:
sage: s = u.cross(v)
sage: s
Vector field u x v on the Euclidean space E^3
sage: s.display()
u x v = (sin(x*y*z)*u_y(x, y, z) - x*u_z(x, y, z)) e_x
+ (-sin(x*y*z)*u_x(x, y, z) - y*u_z(x, y, z)) e_y
+ (x*u_x(x, y, z) + y*u_y(x, y, z)) e_z
We can check the standard formulas expressing the cross product in terms of the components:
sage: 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]])
True
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:
sage: w = E.vector_field(name='w')
sage: w[1] = x*z
sage: w[2] = y*z
sage: w.display()
w = x*z e_x + y*z e_y
The scalar triple product of the vector fields \(u\), \(v\) and \(w\) is obtained as follows:
sage: triple_product = E.scalar_triple_product()
sage: s = triple_product(u, v, w)
sage: s
Scalar field epsilon(u,v,w) on the Euclidean space E^3
sage: s.expr()
-(y*u_x(x, y, z) - x*u_y(x, y, z))*z*sin(x*y*z) - (x^2*u_z(x, y, z)
+ y^2*u_z(x, y, z))*z
Let us check that the scalar triple product of \(u\), \(v\) and \(w\) is \(u\cdot(v\times w)\):
sage: s == u.dot(v.cross(w))
True
How to evaluate the standard differential operators¶
The standard operators \(\mathrm{grad}\), \(\mathrm{div}\), \(\mathrm{curl}\), etc.
involved in vector calculus are accessible as methods on scalar fields and
vector fields (e.g. v.div()). However, to allow for standard mathematical
notations (e.g. div(v)), let us import the functions
grad(), div(),
curl() and
laplacian():
sage: from sage.manifolds.operators import *
Gradient¶
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)\):
sage: F = E.scalar_field(function('f')(x,y,z), name='F')
sage: F.display()
F: E^3 → ℝ
(x, y, z) ↦ f(x, y, z)
The value of \(F\) at a point:
sage: F(p)
f(3, -2, 1)
The gradient of \(F\):
sage: grad(F)
Vector field grad(F) on the Euclidean space E^3
sage: grad(F).display()
grad(F) = d(f)/dx e_x + d(f)/dy e_y + d(f)/dz e_z
sage: norm(grad(F)).display()
|grad(F)|: E^3 → ℝ
(x, y, z) ↦ sqrt((d(f)/dx)^2 + (d(f)/dy)^2 + (d(f)/dz)^2)
Divergence¶
The divergence of the vector field \(u\):
sage: s = div(u)
sage: s.display()
div(u): E^3 → ℝ
(x, y, z) ↦ d(u_x)/dx + d(u_y)/dy + d(u_z)/dz
For \(v\) and \(w\), we have:
sage: div(v).expr()
x*y*cos(x*y*z)
sage: div(w).expr()
2*z
An identity valid for any scalar field \(F\) and any vector field \(u\):
sage: div(F*u) == F*div(u) + u.dot(grad(F))
True
Curl¶
The curl of the vector field \(u\):
sage: s = curl(u)
sage: s
Vector field curl(u) on the Euclidean space E^3
sage: s.display()
curl(u) = (-d(u_y)/dz + d(u_z)/dy) e_x + (d(u_x)/dz - d(u_z)/dx) e_y
+ (-d(u_x)/dy + d(u_y)/dx) e_z
To use the notation rot instead of curl, simply do:
sage: rot = curl
An alternative is:
sage: from sage.manifolds.operators import curl as rot
We have then:
sage: rot(u).display()
curl(u) = (-d(u_y)/dz + d(u_z)/dy) e_x + (d(u_x)/dz - d(u_z)/dx) e_y
+ (-d(u_x)/dy + d(u_y)/dx) e_z
sage: rot(u) == curl(u)
True
For \(v\) and \(w\), we have:
sage: curl(v).display()
curl(v) = x*z*cos(x*y*z) e_x - y*z*cos(x*y*z) e_y + 2 e_z
sage: curl(w).display()
curl(w) = -y e_x + x e_y
The curl of a gradient is always zero:
sage: curl(grad(F)).display()
curl(grad(F)) = 0
The divergence of a curl is always zero:
sage: div(curl(u)).display()
div(curl(u)): E^3 → ℝ
(x, y, z) ↦ 0
An identity valid for any scalar field \(F\) and any vector field \(u\) is
as we can check:
sage: curl(F*u) == grad(F).cross(u) + F*curl(u)
True
Laplacian¶
The Laplacian \(\Delta F\) of a scalar field \(F\) is another scalar field:
sage: s = laplacian(F)
sage: s.display()
Delta(F): E^3 → ℝ
(x, y, z) ↦ d^2(f)/dx^2 + d^2(f)/dy^2 + d^2(f)/dz^2
The following identity holds:
as we can check:
sage: laplacian(F) == div(grad(F))
True
The Laplacian \(\Delta u\) of a vector field \(u\) is another vector field:
sage: Du = laplacian(u)
sage: Du
Vector field Delta(u) on the Euclidean space E^3
whose components are:
sage: Du.display()
Delta(u) = (d^2(u_x)/dx^2 + d^2(u_x)/dy^2 + d^2(u_x)/dz^2) e_x
+ (d^2(u_y)/dx^2 + d^2(u_y)/dy^2 + d^2(u_y)/dz^2) e_y
+ (d^2(u_z)/dx^2 + d^2(u_z)/dy^2 + d^2(u_z)/dz^2) e_z
In the Cartesian frame, the components of \(\Delta u\) are nothing but the (scalar) Laplacians of the components of \(u\), as we can check:
sage: e = E.cartesian_frame()
sage: Du == sum(laplacian(u[[i]])*e[i] for i in E.irange())
True
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:
sage: e[:]
(Vector field e_x on the Euclidean space E^3,
Vector field e_y on the Euclidean space E^3,
Vector field e_z on the Euclidean space E^3)
For the vector fields \(v\) and \(w\), we have:
sage: laplacian(v).display()
Delta(v) = -(x^2*y^2 + (x^2 + y^2)*z^2)*sin(x*y*z) e_z
sage: laplacian(w).display()
Delta(w) = 0
We have:
sage: curl(curl(u)).display()
curl(curl(u)) = (-d^2(u_x)/dy^2 - d^2(u_x)/dz^2 + d^2(u_y)/dxdy
+ d^2(u_z)/dxdz) e_x + (d^2(u_x)/dxdy - d^2(u_y)/dx^2 - d^2(u_y)/dz^2
+ d^2(u_z)/dydz) e_y + (d^2(u_x)/dxdz + d^2(u_y)/dydz - d^2(u_z)/dx^2
- d^2(u_z)/dy^2) e_z
sage: grad(div(u)).display()
grad(div(u)) = (d^2(u_x)/dx^2 + d^2(u_y)/dxdy + d^2(u_z)/dxdz) e_x
+ (d^2(u_x)/dxdy + d^2(u_y)/dy^2 + d^2(u_z)/dydz) e_y
+ (d^2(u_x)/dxdz + d^2(u_y)/dydz + d^2(u_z)/dz^2) e_z
A famous identity is
Let us check it:
sage: curl(curl(u)) == grad(div(u)) - laplacian(u)
True
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)\):
sage: frame = E.cartesian_frame()
sage: frame
Coordinate frame (E^3, (e_x,e_y,e_z))
But this can be changed, thanks to the method
set_name():
sage: frame.set_name('a', indices=('x', 'y', 'z'))
sage: frame
Coordinate frame (E^3, (a_x,a_y,a_z))
sage: v.display()
v = -y a_x + x a_y + sin(x*y*z) a_z
sage: frame.set_name(('hx', 'hy', 'hz'),
....: latex_symbol=(r'\mathbf{\hat{x}}', r'\mathbf{\hat{y}}',
....: r'\mathbf{\hat{z}}'))
sage: frame
Coordinate frame (E^3, (hx,hy,hz))
sage: v.display()
v = -y hx + x hy + sin(x*y*z) hz
Customizing the coordinate symbols¶
The coordinates symbols are defined within the angle brackets <...> at the
construction of the Euclidean space. Above we did:
sage: 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:
sage: E.<X,Y,Z> = EuclideanSpace()
We have then:
sage: E.atlas()
[Chart (E^3, (X, Y, Z))]
sage: E.cartesian_frame()
Coordinate frame (E^3, (e_X,e_Y,e_Z))
sage: v = E.vector_field(-Y, X, sin(X*Y*Z), name='v')
sage: v.display()
v = -Y e_X + X e_Y + sin(X*Y*Z) e_Z
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
(:):
sage: E.<xi,et,ze> = EuclideanSpace(symbols=r"xi:\xi et:\eta ze:\zeta")
sage: E.atlas()
[Chart (E^3, (xi, et, ze))]
sage: v = E.vector_field(-et, xi, sin(xi*et*ze), name='v')
sage: v.display()
v = -et e_xi + xi e_et + sin(et*xi*ze) e_ze