Symbolic and numeric

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

See also

Manual: Elastic moduli

Tutorials: Differentiation and interoperability · Performance of the structured types · Differential operators, symbolically · The elastic sphere, solved symbolically

<!– /crosslinks –>

TensND is generic in its element type. The same code paths run on Float64, on ForwardDiff.Dual, on SymPy's Sym and on Symbolics' Num, which is what lets a derivation be done symbolically and the result evaluated — or differentiated — numerically without rewriting anything.

The four scalar worlds

Element type	Load	Use for
Float64 (and any Real)	using TensND	numerical work
ForwardDiff.Dual	using TensND, ForwardDiff	derivatives and sensitivities
Sym	using TensND, SymPyPythonCall	symbolic derivation, dsolve, solve
Num	using TensND, Symbolics	native Julia CAS, code generation

SymPyPythonCall must be loaded explicitly: it is an optional dependency of TensND, whose SymPy-specific methods (the symbolic simplifications of a Sym, the default coordinates of the coordinate systems) it activates.

Integer and rational coefficients

A structured tensor stores its coefficients in the type its components need. The isotropic and cubic projectors and the Walpole tensors and carry or in their components, and , carry . A fourth-order TensTI built from integers or rationals is therefore stored in floats, and a fourth-order TensISO or a TensCubic built from integers too, while rationals stay exact there; BigInt becomes BigFloat, and every other type is kept as it is:

using TensND
(eltype(TensISO{3}(3, 2)), eltype(TensISO{3}(3 // 1, 2 // 1)),
    eltype(TensTI{4}(1 // 1, 2 // 1, 1 // 1, 1 // 1, 1 // 1, (0, 0, 1))),
    eltype(walpole_basis((0, 0, 1))[3]))

(Float64, Rational{Int64}, Float64, Float64)

A tensor whose components are integers when its coefficients are, a TensTI or a TensISO of order 2 or a TensOrtho, keeps an integer type, and its inverse is in the type of the reciprocal, as inv of an integer matrix is.

using TensND, SymPyPythonCall

k, μ = symbols("k μ", positive = true)
ℂ = TensISO{3}(3k, 2μ)
get_data(inv(ℂ))

(1/(3*k), 1/(2*μ))

The same with Symbolics:

using TensND, Symbolics

Symbolics.@variables kn μn
get_data(inv(TensISO{3}(3kn, 2μn)))

(1 / (3kn), 1 / (2μn))

What is generic, and what is not

	tensor algebra	structured types	projection, fixed axis	projection, optimized	symbolic operators	numerical operators
Float64	✓	✓	✓	✓	—	✓
ForwardDiff.Dual	✓	✓	✓	—	—	✓
Sym	✓	✓	✓	—	✓	—
Num	✓	✓	✓	partial	✓	—

The one genuine restriction is the orientation search: it is a numerical optimization, so it needs a numeric element type. Projection onto a class with a given axis or frame is closed-form and works symbolically.

Symbolic helpers

Six functions apply a symbolic operation elementwise to a tensor, an array or a scalar, and are no-ops on numeric types — so generic code can call them unconditionally:

Function	Does	SymPy	Symbolics
tsimplify	simplify	✓	✓
tfactor	factorize	✓	—
tsubs	substitute	✓	✓
tdiff	differentiate	✓	✓
ttrigsimp	trigonometric simplification	✓	—
texpand_trig	trigonometric expansion	✓	—

tsimplify(3.0), tsimplify(k + k)

(3.0, 2*k)

Being no-ops on Float64 is the point: a function written with tsimplify sprinkled through it runs unchanged, and at full speed, on numeric input.

Differentiating

ForwardDiff composes with the tensor algebra, the structured types and the projections:

using TensND, ForwardDiff, LinearAlgebra

n = [0.0, 0.0, 1.0]
f(α) = get_data(project(get_array(ti_from_constants(10.0, 3.0, 2.5, 12.0, 2.0, [sin(α), 0.0, cos(α)])), :TI, n)[1])[1]

ad = ForwardDiff.derivative(f, 0.3)
fd = (f(0.3 + 1.0e-6) - f(0.3 - 1.0e-6)) / 2.0e-6
(ad, fd)

(-5.32346083364259, -5.323460833928095)

and with the numerical differential operators, in both directions — with respect to the evaluation point and with respect to a parameter carried by the field. See Numerical coordinate systems.

Interoperability with Tensors.jl

TensND stores its data in Tensors.jl containers whenever the shape allows (order 1, 2 or 4, dimension 1 to 3; any other shape keeps a plain array), so the two compose directly and index symmetries are preserved:

using TensND, Tensors, LinearAlgebra

st = SymmetricTensor{2, 3}((i, j) -> Float64(i + j))
t = Tens(st)
typeof(get_array(t))

Tensors.SymmetricTensor{2, 3, Float64, 6}

C4 = SymmetricTensor{4, 3}((i, j, k, l) -> Float64(i + j + k + l))
norm(get_array(from_mandel(mandel(Tens(C4)))) - get_array(Tens(C4)))

2.5121479338940403e-15

Performance notes

Symbolic element types are orders of magnitude slower than numeric ones, and expression size — not operation count — dominates. Three practical consequences:

• simplify early and often. An unsimplified intermediate propagates into every later expression. This is what the rules / tmp_var / to_coords machinery of CoordinateSystem exists for.

• derive symbolically, evaluate numerically. Obtain the closed form once, then substitute with tsubs or generate a numeric function from it.

• prefer a structured type over a dense array whenever the physics allows — the measured gains are one to three orders of magnitude, see Performance of the structured types.

For pointwise numerical work, the pointwise CoordinateSystem avoids the symbolic layer entirely and costs a constant amount per point.