如何在曲线坐标下实现向量微积分

本教程将介绍SageMath在3维欧几里德空间中的一些向量演算功能。相应的工具已经通过 SageManifolds 项目。

The tutorial is also available as a Jupyter notebook, either passive (nbviewer) or interactive (binder).

使用球面坐标

使用球面坐标的步骤 (r,theta,phi) 在欧几里得三维空间内 mathbb{E}^3 ，只要用关键字声明后者就足够了 coordinates='spherical' **

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

多亏了这个符号 <r,th,ph> 在上面的声明中，坐标 (r,theta,phi) 作为三个符号变量立即可用 r ， th 和 ph (不需要通过以下方式申报 var() ，即打字 r, th, ph = var('r th ph') )：

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

此外，坐标LaTeX符号已经设置：：

sage: latex(th)
{\theta}

坐标范围为：

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

mathbb{E}^3 被赋予了 orthonormal 向量标架 (e_r, e_theta, e_phi) 与球面坐标关联：

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

In the above output, (∂/∂r,∂/∂th,∂/∂ph) = left(frac{partial}{partial r}, frac{partial}{partialtheta}, 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. The default frame is the orthonormal one:

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

定义向量场

我们在上定义了一个向量场 mathbb{E}^3 从它在正交化向量坐标系中的分量 (e_r,e_theta,e_phi) **

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

我们可以访问的组件 v 通过方括号运算符：：

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

向量场可以在以下任一点求值 mathbb{E}^3 **

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

我们可以定义一个具有一般分量的矢量场：

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

它在当时的价值 p 然后：：

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

球坐标中的微分算子

标准运算符 mathrm{grad} ， mathrm{div} ， mathrm{curl} 等可作为标量场和矢量场上的方法访问(例如 v.div() )。然而，为了允许使用标准的数学符号(例如 div(v) )，让我们导入函数 grad() ， div() ， curl() 和 laplacian() **

sage: from sage.manifolds.operators import *

梯度

我们首先引入一个标量场，通过它在笛卡尔坐标下的表达式；在这个例子中，我们考虑一个未指定的函数 (r,theta,phi) **

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)

的价值 F 在某一时刻：：

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

的渐变 F **

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

分歧

向量场的散度：：

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)

为 v ，我们有：

sage: div(v).expr()
3

卷曲

向量场的旋度：：

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

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

为 v ，我们有：

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

渐变的旋度始终为零：：

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

卷曲的散度始终为零：：

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

拉普拉斯

标量域的拉普拉斯量：：

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)

向量场的拉普拉斯：：

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

由于此表达式相当长，因此我们可能会逐个请求显示组件：：

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)

我们可以展开每个组件：：

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

作为测试，我们可以检查这些公式是否与维基百科文章中的公式一致 :wikipedia:`Del in cylindrical and spherical coordinates <Del_in_cylindrical_and_spherical_coordinates#Del_formula>` 。

使用柱面坐标

柱面坐标的使用 (rho,phi,z) 在欧几里得空间 mathbb{E}^3 与球面坐标处于相同的基础上。首先，人们只需声明：：

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

然后坐标范围为：：

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

默认向量帧为正交帧 (e_rho,e_phi,e_z) 与柱面坐标相关联：

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

并且可以根据其在该帧中的分量定义矢量场：：

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]

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

柱坐标中的微分算子

sage: from sage.manifolds.operators import *

渐变：：

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

分歧：：

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)

卷曲：：

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

标量域的拉普拉斯量：：

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)

向量场的拉普拉斯：：

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

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

同样，我们可以检查上面的公式是否与维基百科文章中的公式一致 :wikipedia:`Del in cylindrical and spherical coordinates <Del_in_cylindrical_and_spherical_coordinates#Del_formula>` 。

更改坐标

给定一个矢量场在给定坐标系中的表达式，SageMath可以在另一个坐标系中计算其表达式，请参阅教程 如何更改坐标