Differential operators, symbolically

<!– crosslinks: written by .github/scripts/crosslinks.py from docs/crosslinks.toml; do not edit –>

See also

Theory: Curvilinear differential calculus

Manual: Coordinate systems (symbolic) · Coordinate systems (numerical) · Symbolic and numeric · Elastic moduli

Tutorials: Christoffel symbols · Differential operators, numerically · Classical identities of the differential operators · Differentiation and interoperability · Performance of the structured types · The elastic sphere, solved symbolically

<!– /crosslinks –>

The five operators of Curvilinear differential calculus — grad, symgrad, divergence, laplacian, 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 SymPyPythonCall

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, θ = get_coords(Polar)
𝐞ʳ, 𝐞ᶿ = unitvec(Polar)
@set_coorsys Polar

CoordinateSystem{2, Sym{PythonCall.Py}, ExactGeometry{2, Sym{PythonCall.Py}, TensND.TensRotated{1, 2, Sym{PythonCall.Py}, Tensors.Vec{2, Sym{PythonCall.Py}}, RotatedBasis{2, Sym{PythonCall.Py}}}, RotatedBasis{2, Sym{PythonCall.Py}}, OrthogonalBasis{2, Sym{PythonCall.Py}, RotatedBasis{2, Sym{PythonCall.Py}}}}}(ExactGeometry{2, Sym{PythonCall.Py}, TensND.TensRotated{1, 2, Sym{PythonCall.Py}, Tensors.Vec{2, Sym{PythonCall.Py}}, RotatedBasis{2, Sym{PythonCall.Py}}}, RotatedBasis{2, Sym{PythonCall.Py}}, OrthogonalBasis{2, Sym{PythonCall.Py}, RotatedBasis{2, Sym{PythonCall.Py}}}}(Sym{PythonCall.Py}[r, 0], (r, θ), Sym{PythonCall.Py}[cos(θ) -sin(θ); sin(θ) cos(θ)], Sym{PythonCall.Py}[cos(θ) -r*sin(θ); sin(θ) r*cos(θ)], (Sym{PythonCall.Py}[1, 0], Sym{PythonCall.Py}[0, r]), (1, r), (Sym{PythonCall.Py}[1, 0], Sym{PythonCall.Py}[0, 1/r]), (Sym{PythonCall.Py}[1, 0], Sym{PythonCall.Py}[0, 1]), Sym{PythonCall.Py}[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, laplacian reproduces the textbook operator   :

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

  

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

n = symbols("n", integer = true)
simplify(laplacian(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(laplacian(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) - laplacian(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 = get_coords(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 = get_coords(Spherical)
𝐞ᶿ, 𝐞ᵠ, 𝐞ʳˢ = unitvec(Spherical)
@set_coorsys Spherical

get_coords(Spherical), lame_coefficients(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(divergence(𝐞ʳˢ))

(Sym{PythonCall.Py}[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(laplacian(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(divergence(𝛔))

Index placement is a convention

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

Cartesian = coorsys_cartesian()
X = get_coords(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.