How to perform vector calculus in curvilinear coordinates¶

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

Using spherical coordinates¶

To use spherical coordinates $(r, θ, ϕ)$ within the Euclidean 3-space $E^{3}$ , it suffices to declare the latter with the keyword coordinates='spherical':

Sage

sage: E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
sage: E
Euclidean space E^3

Python

>>> from sage.all import *
>>> E = EuclideanSpace(coordinates='spherical', names=('r', 'th', 'ph',)); (r, th, ph,) = E._first_ngens(3)
>>> E
Euclidean space E^3

Sage Live

E.<r,th,ph> = EuclideanSpace(coordinates='spherical')
E

Thanks to the notation <r,th,ph> in the above declaration, the coordinates $(r, θ, ϕ)$ are immediately available as three symbolic variables r, th and ph (there is no need to declare them via var(), i.e. to type r, th, ph = var('r th ph')):

Sage

sage: r is E.spherical_coordinates()[1]
True
sage: (r, th, ph) == E.spherical_coordinates()[:]
True
sage: type(r)
<class 'sage.symbolic.expression.Expression'>

Python

>>> from sage.all import *
>>> r is E.spherical_coordinates()[Integer(1)]
True
>>> (r, th, ph) == E.spherical_coordinates()[:]
True
>>> type(r)
<class 'sage.symbolic.expression.Expression'>

Sage Live

r is E.spherical_coordinates()[1]
(r, th, ph) == E.spherical_coordinates()[:]
type(r)

Moreover, the coordinate LaTeX symbols are already set:

Sage

sage: latex(th)
{\theta}

Python

>>> from sage.all import *
>>> latex(th)
{\theta}

Sage Live

latex(th)

The coordinate ranges are:

Sage

sage: E.spherical_coordinates().coord_range()
r: (0, +oo); th: (0, pi); ph: [0, 2*pi] (periodic)

Python

>>> from sage.all import *
>>> E.spherical_coordinates().coord_range()
r: (0, +oo); th: (0, pi); ph: [0, 2*pi] (periodic)

Sage Live

E.spherical_coordinates().coord_range()

$E^{3}$ is endowed with the orthonormal vector frame $(e_{r}, e_{θ}, e_{ϕ})$ associated with spherical coordinates:

Sage

sage: E.frames()
[Coordinate frame (E^3, (∂/∂r,∂/∂th,∂/∂ph)),
 Vector frame (E^3, (e_r,e_th,e_ph))]

Python

>>> from sage.all import *
>>> E.frames()
[Coordinate frame (E^3, (∂/∂r,∂/∂th,∂/∂ph)),
 Vector frame (E^3, (e_r,e_th,e_ph))]

Sage Live

E.frames()

In the above output, (∂/∂r,∂/∂th,∂/∂ph) = $(\frac{\partial}{\partial r}, \frac{\partial}{\partial θ}, \frac{\partial}{\partial ϕ})$ is the coordinate frame associated with $(r, θ, ϕ)$ ; it is not an orthonormal frame and will not be used below. The default frame is the orthonormal one:

Sage

sage: E.default_frame()
Vector frame (E^3, (e_r,e_th,e_ph))

Python

>>> from sage.all import *
>>> E.default_frame()
Vector frame (E^3, (e_r,e_th,e_ph))

Sage Live

E.default_frame()

Defining a vector field¶

We define a vector field on $E^{3}$ from its components in the orthonormal vector frame $(e_{r}, e_{θ}, e_{ϕ})$ :

Sage

sage: 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')
sage: v.display()
v = (r*sin(2*ph)*sin(th)^2 + r) e_r + r*cos(th)*sin(2*ph)*sin(th) e_th
 + 2*r*cos(ph)^2*sin(th) e_ph

Python

>>> from sage.all import *
>>> v = E.vector_field(r*sin(Integer(2)*ph)*sin(th)**Integer(2) + r,
...                    r*sin(Integer(2)*ph)*sin(th)*cos(th),
...                    Integer(2)*r*cos(ph)**Integer(2)*sin(th), name='v')
>>> v.display()
v = (r*sin(2*ph)*sin(th)^2 + r) e_r + r*cos(th)*sin(2*ph)*sin(th) e_th
 + 2*r*cos(ph)^2*sin(th) e_ph

Sage Live

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:

Sage

sage: v[1]
r*sin(2*ph)*sin(th)^2 + r
sage: v[:]
[r*sin(2*ph)*sin(th)^2 + r, r*cos(th)*sin(2*ph)*sin(th), 2*r*cos(ph)^2*sin(th)]

Python

>>> from sage.all import *
>>> v[Integer(1)]
r*sin(2*ph)*sin(th)^2 + r
>>> v[:]
[r*sin(2*ph)*sin(th)^2 + r, r*cos(th)*sin(2*ph)*sin(th), 2*r*cos(ph)^2*sin(th)]

Sage Live

v[1]
v[:]

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

Sage

sage: p = E((1, pi/2, pi), name='p')
sage: p
Point p on the Euclidean space E^3
sage: p.coordinates()
(1, 1/2*pi, pi)
sage: vp = v.at(p)
sage: vp
Vector v at Point p on the Euclidean space E^3
sage: vp.display()
v = e_r + 2 e_ph

Python

>>> from sage.all import *
>>> p = E((Integer(1), pi/Integer(2), pi), name='p')
>>> p
Point p on the Euclidean space E^3
>>> p.coordinates()
(1, 1/2*pi, pi)
>>> vp = v.at(p)
>>> vp
Vector v at Point p on the Euclidean space E^3
>>> vp.display()
v = e_r + 2 e_ph

Sage Live

p = E((1, pi/2, pi), name='p')
p
p.coordinates()
vp = v.at(p)
vp
vp.display()

We may define a vector field with generic components:

Sage

sage: u = E.vector_field(function('u_r')(r,th,ph),
....:                    function('u_theta')(r,th,ph),
....:                    function('u_phi')(r,th,ph),
....:                    name='u')
sage: u.display()
u = u_r(r, th, ph) e_r + u_theta(r, th, ph) e_th + u_phi(r, th, ph) e_ph
sage: u[:]
[u_r(r, th, ph), u_theta(r, th, ph), u_phi(r, th, ph)]

Python

>>> from sage.all import *
>>> 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 = u_r(r, th, ph) e_r + u_theta(r, th, ph) e_th + u_phi(r, th, ph) e_ph
>>> u[:]
[u_r(r, th, ph), u_theta(r, th, ph), u_phi(r, th, ph)]

Sage Live

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

Its value at the point $p$ is then:

Sage

sage: up = u.at(p)
sage: up.display()
u = u_r(1, 1/2*pi, pi) e_r + u_theta(1, 1/2*pi, pi) e_th
 + u_phi(1, 1/2*pi, pi) e_ph

Python

>>> from sage.all import *
>>> up = u.at(p)
>>> up.display()
u = u_r(1, 1/2*pi, pi) e_r + u_theta(1, 1/2*pi, pi) e_th
 + u_phi(1, 1/2*pi, pi) e_ph

Sage Live

up = u.at(p)
up.display()

Differential operators in spherical coordinates¶

The standard operators $grad$ , $div$ , $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

sage: from sage.manifolds.operators import *

Python

>>> from sage.all import *
>>> from sage.manifolds.operators import *

Sage Live

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 a unspecified function of $(r, θ, ϕ)$ :

Sage

sage: F = E.scalar_field(function('f')(r,th,ph), name='F')
sage: F.display()
F: E^3 → ℝ
   (r, th, ph) ↦ f(r, th, ph)

Python

>>> from sage.all import *
>>> F = E.scalar_field(function('f')(r,th,ph), name='F')
>>> F.display()
F: E^3 → ℝ
   (r, th, ph) ↦ f(r, th, ph)

Sage Live

F = E.scalar_field(function('f')(r,th,ph), name='F')
F.display()

The value of $F$ at a point:

Sage

sage: F(p)
f(1, 1/2*pi, pi)

Python

>>> from sage.all import *
>>> F(p)
f(1, 1/2*pi, pi)

Sage Live

F(p)

The gradient of $F$ :

Sage

sage: grad(F)
Vector field grad(F) on the Euclidean space E^3
sage: grad(F).display()
grad(F) = d(f)/dr e_r + d(f)/dth/r e_th + d(f)/dph/(r*sin(th)) e_ph
sage: norm(grad(F)).display()
|grad(F)|: E^3 → ℝ
   (r, th, ph) ↦ sqrt((r^2*(d(f)/dr)^2 + (d(f)/dth)^2)*sin(th)^2
    + (d(f)/dph)^2)/(r*sin(th))

Python

>>> from sage.all import *
>>> grad(F)
Vector field grad(F) on the Euclidean space E^3
>>> grad(F).display()
grad(F) = d(f)/dr e_r + d(f)/dth/r e_th + d(f)/dph/(r*sin(th)) e_ph
>>> norm(grad(F)).display()
|grad(F)|: E^3 → ℝ
   (r, th, ph) ↦ sqrt((r^2*(d(f)/dr)^2 + (d(f)/dth)^2)*sin(th)^2
    + (d(f)/dph)^2)/(r*sin(th))

Sage Live

grad(F)
grad(F).display()
norm(grad(F)).display()

Divergence¶

The divergence of a vector field:

Sage

sage: s = div(u)
sage: s.display()
div(u): E^3 → ℝ
   (r, th, ph) ↦ ((r*d(u_r)/dr + 2*u_r(r, th, ph)
    + d(u_theta)/dth)*sin(th) + cos(th)*u_theta(r, th, ph)
    + d(u_phi)/dph)/(r*sin(th))
sage: s.expr().expand()
2*u_r(r, th, ph)/r + cos(th)*u_theta(r, th, ph)/(r*sin(th))
 + diff(u_theta(r, th, ph), th)/r + diff(u_phi(r, th, ph), ph)/(r*sin(th))
 + diff(u_r(r, th, ph), r)

Python

>>> from sage.all import *
>>> s = div(u)
>>> s.display()
div(u): E^3 → ℝ
   (r, th, ph) ↦ ((r*d(u_r)/dr + 2*u_r(r, th, ph)
    + d(u_theta)/dth)*sin(th) + cos(th)*u_theta(r, th, ph)
    + d(u_phi)/dph)/(r*sin(th))
>>> s.expr().expand()
2*u_r(r, th, ph)/r + cos(th)*u_theta(r, th, ph)/(r*sin(th))
 + diff(u_theta(r, th, ph), th)/r + diff(u_phi(r, th, ph), ph)/(r*sin(th))
 + diff(u_r(r, th, ph), r)

Sage Live

s = div(u)
s.display()
s.expr().expand()

For $v$ , we have:

Sage

sage: div(v).expr()
3

Python

>>> from sage.all import *
>>> div(v).expr()
3

Sage Live

div(v).expr()

Curl¶

The curl of a vector field:

Sage

sage: s = curl(u)
sage: s
Vector field curl(u) on the Euclidean space E^3

Python

>>> from sage.all import *
>>> s = curl(u)
>>> s
Vector field curl(u) on the Euclidean space E^3

Sage Live

s = curl(u)
s

Sage

sage: s.display()
curl(u) = (cos(th)*u_phi(r, th, ph) + sin(th)*d(u_phi)/dth
 - d(u_theta)/dph)/(r*sin(th)) e_r - ((r*d(u_phi)/dr + u_phi(r, th, ph))*sin(th)
 - d(u_r)/dph)/(r*sin(th)) e_th + (r*d(u_theta)/dr + u_theta(r, th, ph)
 - d(u_r)/dth)/r e_ph

Python

>>> from sage.all import *
>>> s.display()
curl(u) = (cos(th)*u_phi(r, th, ph) + sin(th)*d(u_phi)/dth
 - d(u_theta)/dph)/(r*sin(th)) e_r - ((r*d(u_phi)/dr + u_phi(r, th, ph))*sin(th)
 - d(u_r)/dph)/(r*sin(th)) e_th + (r*d(u_theta)/dr + u_theta(r, th, ph)
 - d(u_r)/dth)/r e_ph

Sage Live

s.display()

For $v$ , we have:

Sage

sage: curl(v).display()
curl(v) = 2*cos(th) e_r - 2*sin(th) e_th

Python

>>> from sage.all import *
>>> curl(v).display()
curl(v) = 2*cos(th) e_r - 2*sin(th) e_th

Sage Live

curl(v).display()

The curl of a gradient is always zero:

Sage

sage: curl(grad(F)).display()
curl(grad(F)) = 0

Python

>>> from sage.all import *
>>> curl(grad(F)).display()
curl(grad(F)) = 0

Sage Live

curl(grad(F)).display()

The divergence of a curl is always zero:

Sage

sage: div(curl(u)).display()
div(curl(u)): E^3 → ℝ
   (r, th, ph) ↦ 0

Python

>>> from sage.all import *
>>> div(curl(u)).display()
div(curl(u)): E^3 → ℝ
   (r, th, ph) ↦ 0

Sage Live

div(curl(u)).display()

Laplacian¶

The Laplacian of a scalar field:

Sage

sage: s = laplacian(F)
sage: s.display()
Delta(F): E^3 → ℝ
   (r, th, ph) ↦ ((r^2*d^2(f)/dr^2 + 2*r*d(f)/dr
    + d^2(f)/dth^2)*sin(th)^2 + cos(th)*sin(th)*d(f)/dth
    + d^2(f)/dph^2)/(r^2*sin(th)^2)
sage: s.expr().expand()
2*diff(f(r, th, ph), r)/r + cos(th)*diff(f(r, th, ph), th)/(r^2*sin(th))
 + diff(f(r, th, ph), th, th)/r^2 + diff(f(r, th, ph), ph, ph)/(r^2*sin(th)^2)
 + diff(f(r, th, ph), r, r)

Python

>>> from sage.all import *
>>> s = laplacian(F)
>>> s.display()
Delta(F): E^3 → ℝ
   (r, th, ph) ↦ ((r^2*d^2(f)/dr^2 + 2*r*d(f)/dr
    + d^2(f)/dth^2)*sin(th)^2 + cos(th)*sin(th)*d(f)/dth
    + d^2(f)/dph^2)/(r^2*sin(th)^2)
>>> s.expr().expand()
2*diff(f(r, th, ph), r)/r + cos(th)*diff(f(r, th, ph), th)/(r^2*sin(th))
 + diff(f(r, th, ph), th, th)/r^2 + diff(f(r, th, ph), ph, ph)/(r^2*sin(th)^2)
 + diff(f(r, th, ph), r, r)

Sage Live

s = laplacian(F)
s.display()
s.expr().expand()

The Laplacian of a vector field:

Sage

sage: Du = laplacian(u)
sage: Du.display()
Delta(u) = ((r^2*d^2(u_r)/dr^2 + 2*r*d(u_r)/dr - 2*u_r(r, th, ph)
 + d^2(u_r)/dth^2 - 2*d(u_theta)/dth)*sin(th)^2 - ((2*u_theta(r, th, ph)
 - d(u_r)/dth)*cos(th) + 2*d(u_phi)/dph)*sin(th) + d^2(u_r)/dph^2)/(r^2*sin(th)^2) e_r
 + ((r^2*d^2(u_theta)/dr^2 + 2*r*d(u_theta)/dr + 2*d(u_r)/dth + d^2(u_theta)/dth^2)*sin(th)^2
 + cos(th)*sin(th)*d(u_theta)/dth - 2*cos(th)*d(u_phi)/dph - u_theta(r, th, ph)
 + d^2(u_theta)/dph^2)/(r^2*sin(th)^2) e_th
 + ((r^2*d^2(u_phi)/dr^2 + 2*r*d(u_phi)/dr
 + d^2(u_phi)/dth^2)*sin(th)^2 + (cos(th)*d(u_phi)/dth + 2*d(u_r)/dph)*sin(th)
 + 2*cos(th)*d(u_theta)/dph - u_phi(r, th, ph) + d^2(u_phi)/dph^2)/(r^2*sin(th)^2) e_ph

Python

>>> from sage.all import *
>>> Du = laplacian(u)
>>> Du.display()
Delta(u) = ((r^2*d^2(u_r)/dr^2 + 2*r*d(u_r)/dr - 2*u_r(r, th, ph)
 + d^2(u_r)/dth^2 - 2*d(u_theta)/dth)*sin(th)^2 - ((2*u_theta(r, th, ph)
 - d(u_r)/dth)*cos(th) + 2*d(u_phi)/dph)*sin(th) + d^2(u_r)/dph^2)/(r^2*sin(th)^2) e_r
 + ((r^2*d^2(u_theta)/dr^2 + 2*r*d(u_theta)/dr + 2*d(u_r)/dth + d^2(u_theta)/dth^2)*sin(th)^2
 + cos(th)*sin(th)*d(u_theta)/dth - 2*cos(th)*d(u_phi)/dph - u_theta(r, th, ph)
 + d^2(u_theta)/dph^2)/(r^2*sin(th)^2) e_th
 + ((r^2*d^2(u_phi)/dr^2 + 2*r*d(u_phi)/dr
 + d^2(u_phi)/dth^2)*sin(th)^2 + (cos(th)*d(u_phi)/dth + 2*d(u_r)/dph)*sin(th)
 + 2*cos(th)*d(u_theta)/dph - u_phi(r, th, ph) + d^2(u_phi)/dph^2)/(r^2*sin(th)^2) e_ph

Sage Live

Du = laplacian(u)
Du.display()

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

Sage

sage: Du.display_comp()
Delta(u)^1 = ((r^2*d^2(u_r)/dr^2 + 2*r*d(u_r)/dr - 2*u_r(r, th, ph) + d^2(u_r)/dth^2
 - 2*d(u_theta)/dth)*sin(th)^2 - ((2*u_theta(r, th, ph) - d(u_r)/dth)*cos(th)
 + 2*d(u_phi)/dph)*sin(th) + d^2(u_r)/dph^2)/(r^2*sin(th)^2)
Delta(u)^2 = ((r^2*d^2(u_theta)/dr^2 + 2*r*d(u_theta)/dr + 2*d(u_r)/dth
 + d^2(u_theta)/dth^2)*sin(th)^2 + cos(th)*sin(th)*d(u_theta)/dth
 - 2*cos(th)*d(u_phi)/dph - u_theta(r, th, ph) + d^2(u_theta)/dph^2)/(r^2*sin(th)^2)
Delta(u)^3 = ((r^2*d^2(u_phi)/dr^2 + 2*r*d(u_phi)/dr + d^2(u_phi)/dth^2)*sin(th)^2
 + (cos(th)*d(u_phi)/dth + 2*d(u_r)/dph)*sin(th) + 2*cos(th)*d(u_theta)/dph
 - u_phi(r, th, ph) + d^2(u_phi)/dph^2)/(r^2*sin(th)^2)

Python

>>> from sage.all import *
>>> Du.display_comp()
Delta(u)^1 = ((r^2*d^2(u_r)/dr^2 + 2*r*d(u_r)/dr - 2*u_r(r, th, ph) + d^2(u_r)/dth^2
 - 2*d(u_theta)/dth)*sin(th)^2 - ((2*u_theta(r, th, ph) - d(u_r)/dth)*cos(th)
 + 2*d(u_phi)/dph)*sin(th) + d^2(u_r)/dph^2)/(r^2*sin(th)^2)
Delta(u)^2 = ((r^2*d^2(u_theta)/dr^2 + 2*r*d(u_theta)/dr + 2*d(u_r)/dth
 + d^2(u_theta)/dth^2)*sin(th)^2 + cos(th)*sin(th)*d(u_theta)/dth
 - 2*cos(th)*d(u_phi)/dph - u_theta(r, th, ph) + d^2(u_theta)/dph^2)/(r^2*sin(th)^2)
Delta(u)^3 = ((r^2*d^2(u_phi)/dr^2 + 2*r*d(u_phi)/dr + d^2(u_phi)/dth^2)*sin(th)^2
 + (cos(th)*d(u_phi)/dth + 2*d(u_r)/dph)*sin(th) + 2*cos(th)*d(u_theta)/dph
 - u_phi(r, th, ph) + d^2(u_phi)/dph^2)/(r^2*sin(th)^2)

Sage Live

Du.display_comp()

We may expand each component:

Sage

sage: for i in E.irange():
....:     s = Du[i].expand()
sage: Du.display_comp()
Delta(u)^1 = 2*d(u_r)/dr/r - 2*u_r(r, th, ph)/r^2
 - 2*cos(th)*u_theta(r, th, ph)/(r^2*sin(th)) + cos(th)*d(u_r)/dth/(r^2*sin(th))
 + d^2(u_r)/dth^2/r^2 - 2*d(u_theta)/dth/r^2 - 2*d(u_phi)/dph/(r^2*sin(th))
 + d^2(u_r)/dph^2/(r^2*sin(th)^2) + d^2(u_r)/dr^2
Delta(u)^2 = 2*d(u_theta)/dr/r + 2*d(u_r)/dth/r^2 + cos(th)*d(u_theta)/dth/(r^2*sin(th))
 + d^2(u_theta)/dth^2/r^2 - 2*cos(th)*d(u_phi)/dph/(r^2*sin(th)^2)
 - u_theta(r, th, ph)/(r^2*sin(th)^2) + d^2(u_theta)/dph^2/(r^2*sin(th)^2)
 + d^2(u_theta)/dr^2
Delta(u)^3 = 2*d(u_phi)/dr/r + cos(th)*d(u_phi)/dth/(r^2*sin(th))
 + d^2(u_phi)/dth^2/r^2 + 2*d(u_r)/dph/(r^2*sin(th))
 + 2*cos(th)*d(u_theta)/dph/(r^2*sin(th)^2) - u_phi(r, th, ph)/(r^2*sin(th)^2)
 + d^2(u_phi)/dph^2/(r^2*sin(th)^2) + d^2(u_phi)/dr^2

Python

>>> from sage.all import *
>>> for i in E.irange():
...     s = Du[i].expand()
>>> Du.display_comp()
Delta(u)^1 = 2*d(u_r)/dr/r - 2*u_r(r, th, ph)/r^2
 - 2*cos(th)*u_theta(r, th, ph)/(r^2*sin(th)) + cos(th)*d(u_r)/dth/(r^2*sin(th))
 + d^2(u_r)/dth^2/r^2 - 2*d(u_theta)/dth/r^2 - 2*d(u_phi)/dph/(r^2*sin(th))
 + d^2(u_r)/dph^2/(r^2*sin(th)^2) + d^2(u_r)/dr^2
Delta(u)^2 = 2*d(u_theta)/dr/r + 2*d(u_r)/dth/r^2 + cos(th)*d(u_theta)/dth/(r^2*sin(th))
 + d^2(u_theta)/dth^2/r^2 - 2*cos(th)*d(u_phi)/dph/(r^2*sin(th)^2)
 - u_theta(r, th, ph)/(r^2*sin(th)^2) + d^2(u_theta)/dph^2/(r^2*sin(th)^2)
 + d^2(u_theta)/dr^2
Delta(u)^3 = 2*d(u_phi)/dr/r + cos(th)*d(u_phi)/dth/(r^2*sin(th))
 + d^2(u_phi)/dth^2/r^2 + 2*d(u_r)/dph/(r^2*sin(th))
 + 2*cos(th)*d(u_theta)/dph/(r^2*sin(th)^2) - u_phi(r, th, ph)/(r^2*sin(th)^2)
 + d^2(u_phi)/dph^2/(r^2*sin(th)^2) + d^2(u_phi)/dr^2

Sage Live

for i in E.irange():
    s = Du[i].expand()
Du.display_comp()

As a test, we may check that these formulas coincide with those of Wikipedia’s article Del in cylindrical and spherical coordinates.

Using cylindrical coordinates¶

The use of cylindrical coordinates $(ρ, ϕ, z)$ in the Euclidean space $E^{3}$ is on the same footing as that of spherical coordinates. To start with, one has simply to declare:

Sage

sage: E.<rh,ph,z> = EuclideanSpace(coordinates='cylindrical')

Python

>>> from sage.all import *
>>> E = EuclideanSpace(coordinates='cylindrical', names=('rh', 'ph', 'z',)); (rh, ph, z,) = E._first_ngens(3)

Sage Live

E.<rh,ph,z> = EuclideanSpace(coordinates='cylindrical')

The coordinate ranges are then:

Sage

sage: E.cylindrical_coordinates().coord_range()
rh: (0, +oo); ph: [0, 2*pi] (periodic); z: (-oo, +oo)

Python

>>> from sage.all import *
>>> E.cylindrical_coordinates().coord_range()
rh: (0, +oo); ph: [0, 2*pi] (periodic); z: (-oo, +oo)

Sage Live

E.cylindrical_coordinates().coord_range()

The default vector frame is the orthonormal frame $(e_{ρ}, e_{ϕ}, e_{z})$ associated with cylindrical coordinates:

Sage

sage: E.default_frame()
Vector frame (E^3, (e_rh,e_ph,e_z))

Python

>>> from sage.all import *
>>> E.default_frame()
Vector frame (E^3, (e_rh,e_ph,e_z))

Sage Live

E.default_frame()

and one may define vector fields from their components in that frame:

Sage

sage: v = E.vector_field(rh*(1+sin(2*ph)), 2*rh*cos(ph)^2, z,
....:                    name='v')
sage: v.display()
v = rh*(sin(2*ph) + 1) e_rh + 2*rh*cos(ph)^2 e_ph + z e_z
sage: v[:]
[rh*(sin(2*ph) + 1), 2*rh*cos(ph)^2, z]

Python

>>> from sage.all import *
>>> v = E.vector_field(rh*(Integer(1)+sin(Integer(2)*ph)), Integer(2)*rh*cos(ph)**Integer(2), z,
...                    name='v')
>>> v.display()
v = rh*(sin(2*ph) + 1) e_rh + 2*rh*cos(ph)^2 e_ph + z e_z
>>> v[:]
[rh*(sin(2*ph) + 1), 2*rh*cos(ph)^2, z]

Sage Live

v = E.vector_field(rh*(1+sin(2*ph)), 2*rh*cos(ph)^2, z,
                   name='v')
v.display()
v[:]

Sage

sage: u = E.vector_field(function('u_rho')(rh,ph,z),
....:                    function('u_phi')(rh,ph,z),
....:                    function('u_z')(rh,ph,z),
....:                    name='u')
sage: u.display()
u = u_rho(rh, ph, z) e_rh + u_phi(rh, ph, z) e_ph + u_z(rh, ph, z) e_z
sage: u[:]
[u_rho(rh, ph, z), u_phi(rh, ph, z), u_z(rh, ph, z)]

Python

>>> from sage.all import *
>>> 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 = u_rho(rh, ph, z) e_rh + u_phi(rh, ph, z) e_ph + u_z(rh, ph, z) e_z
>>> u[:]
[u_rho(rh, ph, z), u_phi(rh, ph, z), u_z(rh, ph, z)]

Sage Live

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

Differential operators in cylindrical coordinates¶

Sage

sage: from sage.manifolds.operators import *

Python

>>> from sage.all import *
>>> from sage.manifolds.operators import *

Sage Live

from sage.manifolds.operators import *

The gradient:

Sage

sage: F = E.scalar_field(function('f')(rh,ph,z), name='F')
sage: F.display()
F: E^3 → ℝ
   (rh, ph, z) ↦ f(rh, ph, z)
sage: grad(F)
Vector field grad(F) on the Euclidean space E^3
sage: grad(F).display()
grad(F) = d(f)/drh e_rh + d(f)/dph/rh e_ph + d(f)/dz e_z

Python

>>> from sage.all import *
>>> F = E.scalar_field(function('f')(rh,ph,z), name='F')
>>> F.display()
F: E^3 → ℝ
   (rh, ph, z) ↦ f(rh, ph, z)
>>> grad(F)
Vector field grad(F) on the Euclidean space E^3
>>> grad(F).display()
grad(F) = d(f)/drh e_rh + d(f)/dph/rh e_ph + d(f)/dz e_z

Sage Live

F = E.scalar_field(function('f')(rh,ph,z), name='F')
F.display()
grad(F)
grad(F).display()

The divergence:

Sage

sage: s = div(u)
sage: s.display()
div(u): E^3 → ℝ
   (rh, ph, z) ↦ (rh*d(u_rho)/drh + rh*d(u_z)/dz + u_rho(rh, ph, z) + d(u_phi)/dph)/rh
sage: s.expr().expand()
u_rho(rh, ph, z)/rh + diff(u_phi(rh, ph, z), ph)/rh + diff(u_rho(rh, ph, z), rh)
 + diff(u_z(rh, ph, z), z)

Python

>>> from sage.all import *
>>> s = div(u)
>>> s.display()
div(u): E^3 → ℝ
   (rh, ph, z) ↦ (rh*d(u_rho)/drh + rh*d(u_z)/dz + u_rho(rh, ph, z) + d(u_phi)/dph)/rh
>>> s.expr().expand()
u_rho(rh, ph, z)/rh + diff(u_phi(rh, ph, z), ph)/rh + diff(u_rho(rh, ph, z), rh)
 + diff(u_z(rh, ph, z), z)

Sage Live

s = div(u)
s.display()
s.expr().expand()

The curl:

Sage

sage: s = curl(u)
sage: s
Vector field curl(u) on the Euclidean space E^3
sage: s.display()
curl(u) = -(rh*d(u_phi)/dz - d(u_z)/dph)/rh e_rh + (d(u_rho)/dz - d(u_z)/drh) e_ph
 + (rh*d(u_phi)/drh + u_phi(rh, ph, z) - d(u_rho)/dph)/rh e_z

Python

>>> from sage.all import *
>>> s = curl(u)
>>> s
Vector field curl(u) on the Euclidean space E^3
>>> s.display()
curl(u) = -(rh*d(u_phi)/dz - d(u_z)/dph)/rh e_rh + (d(u_rho)/dz - d(u_z)/drh) e_ph
 + (rh*d(u_phi)/drh + u_phi(rh, ph, z) - d(u_rho)/dph)/rh e_z

Sage Live

s = curl(u)
s
s.display()

The Laplacian of a scalar field:

Sage

sage: s = laplacian(F)
sage: s.display()
Delta(F): E^3 → ℝ
   (rh, ph, z) ↦ (rh^2*d^2(f)/drh^2 + rh^2*d^2(f)/dz^2 + rh*d(f)/drh
    + d^2(f)/dph^2)/rh^2
sage: s.expr().expand()
diff(f(rh, ph, z), rh)/rh + diff(f(rh, ph, z), ph, ph)/rh^2
 + diff(f(rh, ph, z), rh, rh) + diff(f(rh, ph, z), z, z)

Python

>>> from sage.all import *
>>> s = laplacian(F)
>>> s.display()
Delta(F): E^3 → ℝ
   (rh, ph, z) ↦ (rh^2*d^2(f)/drh^2 + rh^2*d^2(f)/dz^2 + rh*d(f)/drh
    + d^2(f)/dph^2)/rh^2
>>> s.expr().expand()
diff(f(rh, ph, z), rh)/rh + diff(f(rh, ph, z), ph, ph)/rh^2
 + diff(f(rh, ph, z), rh, rh) + diff(f(rh, ph, z), z, z)

Sage Live

s = laplacian(F)
s.display()
s.expr().expand()

The Laplacian of a vector field:

Sage

sage: Du = laplacian(u)
sage: Du.display()
Delta(u) = (rh^2*d^2(u_rho)/drh^2 + rh^2*d^2(u_rho)/dz^2 + rh*d(u_rho)/drh
 - u_rho(rh, ph, z) - 2*d(u_phi)/dph + d^2(u_rho)/dph^2)/rh^2 e_rh
 + (rh^2*d^2(u_phi)/drh^2 + rh^2*d^2(u_phi)/dz^2 + rh*d(u_phi)/drh
 - u_phi(rh, ph, z) + d^2(u_phi)/dph^2 + 2*d(u_rho)/dph)/rh^2 e_ph
 + (rh^2*d^2(u_z)/drh^2 + rh^2*d^2(u_z)/dz^2 + rh*d(u_z)/drh
 + d^2(u_z)/dph^2)/rh^2 e_z

Python

>>> from sage.all import *
>>> Du = laplacian(u)
>>> Du.display()
Delta(u) = (rh^2*d^2(u_rho)/drh^2 + rh^2*d^2(u_rho)/dz^2 + rh*d(u_rho)/drh
 - u_rho(rh, ph, z) - 2*d(u_phi)/dph + d^2(u_rho)/dph^2)/rh^2 e_rh
 + (rh^2*d^2(u_phi)/drh^2 + rh^2*d^2(u_phi)/dz^2 + rh*d(u_phi)/drh
 - u_phi(rh, ph, z) + d^2(u_phi)/dph^2 + 2*d(u_rho)/dph)/rh^2 e_ph
 + (rh^2*d^2(u_z)/drh^2 + rh^2*d^2(u_z)/dz^2 + rh*d(u_z)/drh
 + d^2(u_z)/dph^2)/rh^2 e_z

Sage Live

Du = laplacian(u)
Du.display()

Sage

sage: for i in E.irange():
....:     s = Du[i].expand()
sage: Du.display_comp()
Delta(u)^1 = d(u_rho)/drh/rh - u_rho(rh, ph, z)/rh^2 - 2*d(u_phi)/dph/rh^2
 + d^2(u_rho)/dph^2/rh^2 + d^2(u_rho)/drh^2 + d^2(u_rho)/dz^2
Delta(u)^2 = d(u_phi)/drh/rh - u_phi(rh, ph, z)/rh^2 + d^2(u_phi)/dph^2/rh^2
 + 2*d(u_rho)/dph/rh^2 + d^2(u_phi)/drh^2 + d^2(u_phi)/dz^2
Delta(u)^3 = d(u_z)/drh/rh + d^2(u_z)/dph^2/rh^2 + d^2(u_z)/drh^2 + d^2(u_z)/dz^2

Python

>>> from sage.all import *
>>> for i in E.irange():
...     s = Du[i].expand()
>>> Du.display_comp()
Delta(u)^1 = d(u_rho)/drh/rh - u_rho(rh, ph, z)/rh^2 - 2*d(u_phi)/dph/rh^2
 + d^2(u_rho)/dph^2/rh^2 + d^2(u_rho)/drh^2 + d^2(u_rho)/dz^2
Delta(u)^2 = d(u_phi)/drh/rh - u_phi(rh, ph, z)/rh^2 + d^2(u_phi)/dph^2/rh^2
 + 2*d(u_rho)/dph/rh^2 + d^2(u_phi)/drh^2 + d^2(u_phi)/dz^2
Delta(u)^3 = d(u_z)/drh/rh + d^2(u_z)/dph^2/rh^2 + d^2(u_z)/drh^2 + d^2(u_z)/dz^2

Sage Live

for i in E.irange():
    s = Du[i].expand()
Du.display_comp()

Again, we may check that the above formulas coincide with those of Wikipedia’s article Del in cylindrical and spherical coordinates.

Changing coordinates¶

Given the expression of a vector field in a given coordinate system, SageMath can compute its expression in another coordinate system, see the tutorial How to change coordinates