Differential operators, symbolically

The five operators of Curvilinear differential calculus — GRAD, SYMGRAD, DIV, LAPLACE, HESS — applied in polar, cylindrical and spherical coordinates, with exact SymPy derivatives. The results are closed-form expressions, valid everywhere, not values at a point.

The definitions, once:

using TensND
using LinearAlgebra
using SymPy

Setting a default coordinate system

@set_coorsys makes a system the default, so the operators can be called with one argument.

Polar = coorsys_polar()
r, θ = getcoords(Polar)
𝐞ʳ, 𝐞ᶿ = unitvec(Polar)
@set_coorsys Polar

CoorSystemSym{2, Sym{PyCall.PyObject}, TensND.TensRotated{1, 2, Sym{PyCall.PyObject}, Tensors.Vec{2, Sym{PyCall.PyObject}}}, RotatedBasis{2, Sym{PyCall.PyObject}}, OrthogonalBasis{2, Sym{PyCall.PyObject}}}(TensND.ChartCore{2, Sym{PyCall.PyObject}, TensND.TensRotated{1, 2, Sym{PyCall.PyObject}, Tensors.Vec{2, Sym{PyCall.PyObject}}}, RotatedBasis{2, Sym{PyCall.PyObject}}, OrthogonalBasis{2, Sym{PyCall.PyObject}}}(Sym{PyCall.PyObject}[r, 0], (r, θ), Sym{PyCall.PyObject}[cos(θ) -sin(θ); sin(θ) cos(θ)], Sym{PyCall.PyObject}[cos(θ) -r*sin(θ); sin(θ) r*cos(θ)], (Sym{PyCall.PyObject}[1, 0], Sym{PyCall.PyObject}[0, r]), (1, r), (Sym{PyCall.PyObject}[1, 0], Sym{PyCall.PyObject}[0, 1/r]), (Sym{PyCall.PyObject}[1, 0], Sym{PyCall.PyObject}[0, 1]), Sym{PyCall.PyObject}[0 0; 0 -r;;; 0 1/r; 1/r 0], (), (), Dict{Any, Any}(), Dict{Any, Any}(), Dict{Any, Any}()))

The polar Laplacian

Applied to an arbitrary function, LAPLACE reproduces the textbook operator   :

f = SymFunction("f", real = true)(r, θ)
LAPLACE(f)

  

Harmonic functions in the plane are and ; their Laplacian vanishes identically for symbolic :

n = symbols("n", integer = true)
simplify(LAPLACE(r^n * cos(n * θ)))

while a mismatched pair of exponents does not vanish, and the residual is the familiar  :

m = symbols("m", integer = true)
simplify(LAPLACE(r^n * sin(m * θ)) / (r^(n - 2) * sin(m * θ)))

 

The Hessian

  , and its trace is the Laplacian:

H = simplify(HESS(r^n))

simplify(tr(H) - LAPLACE(r^n))

Cylindrical coordinates: an axisymmetric strain tensor

SYMGRAD applied to a displacement field is directly the linearized strain tensor. For an axisymmetric field   , the hoop strain   appears from the Christoffel terms alone — there is no derivative anywhere in the field.

Cylindrical = coorsys_cylindrical()
rc, θc, zc = getcoords(Cylindrical)
𝐞ʳᶜ, 𝐞ᶿᶜ, 𝐞ᶻᶜ = unitvec(Cylindrical)
@set_coorsys Cylindrical

ξʳ = SymFunction("ξʳ", real = true)(rc, zc)
ξᶻ = SymFunction("ξᶻ", real = true)(rc, zc)
𝛏 = ξʳ * 𝐞ʳᶜ + ξᶻ * 𝐞ᶻᶜ

𝛆 = tsimplify(SYMGRAD(𝛏))

Read off the hoop component:

get_array(𝛆)[2, 2]

Spherical coordinates

Recall the unusual ordering , chosen so that    gives the canonical basis in the canonical order.

Spherical = coorsys_spherical()
θs, ϕs, rs = getcoords(Spherical)
𝐞ᶿ, 𝐞ᵠ, 𝐞ʳˢ = unitvec(Spherical)
@set_coorsys Spherical

getcoords(Spherical), Lame(Spherical)

((θ, ϕ, r), (r, r*sin(θ), 1))

The gradient of the radius is the radial unit vector, and the divergence of that vector is :

tsimplify(GRAD(rs)), tsimplify(DIV(𝐞ʳˢ))

(Sym{PyCall.PyObject}[0, 0, 1], 2/r)

The gradient of is the transverse projector divided by — the relation that makes a sphere's curvature , see Submanifolds:

tsimplify(GRAD(𝐞ʳˢ) - (𝐞ᶿ ⊗ 𝐞ᶿ + 𝐞ᵠ ⊗ 𝐞ᵠ) / rs)

is the kernel of the three-dimensional Laplace equation:

tsimplify(LAPLACE(1 / rs))

Radial equilibrium of a spherical stress state

A diagonal stress field in the spherical frame. Its divergence gives the equilibrium equation of a spherically symmetric problem — the   term coming entirely from the connection.

σʳʳ = SymFunction("σʳʳ", real = true)(rs)
σᶿᶿ = SymFunction("σᶿᶿ", real = true)(rs)
σᵠᵠ = SymFunction("σᵠᵠ", real = true)(rs)
𝛔 = σʳʳ * 𝐞ʳˢ ⊗ 𝐞ʳˢ + σᶿᶿ * 𝐞ᶿ ⊗ 𝐞ᶿ + σᵠᵠ * 𝐞ᵠ ⊗ 𝐞ᵠ

simplify(DIV(𝛔))

Index placement is a convention

GRAD appends the derivative index on the right, so  , and DIV contracts the last index,   . In Cartesian coordinates, where all Christoffel symbols vanish, this is easy to read off:

Cartesian = coorsys_cartesian()
X = getcoords(Cartesian)
E = unitvec(Cartesian)
@set_coorsys Cartesian

v = sum(SymFunction("v$i", real = true)(X...) * E[i] for i in 1:3)
Gv = get_array(GRAD(v))

(Gv[1, 2], diff(SymFunction("v1", real = true)(X...), X[2]))

(Derivative(v1(x, y, z), y), Derivative(v1(x, y, z), y))

The entry is , not . A library using the opposite convention differs by a transpose.

---

This page was generated using Literate.jl.