Calculus

import it.grypho.scala.leonardo.core.*
import it.grypho.scala.leonardo.scalar.*

val x = _Variable("x")
val env = new Environment()

Symbolic differentiation

derive(e, v) applies the standard differentiation rules and returns a symbolic expression. The chain rule, product rule, and quotient rule are all handled:

derive(Power(x, _Number(3.0)), x).toString
// res0: String = "(3.0 * (x ^ 2.0))"

// d/dx sin(x^2) = cos(x^2) * 2x   (chain rule)
derive(Sin(Power(x, _Number(2.0))), x).toString
// res1: String = "(cos((x ^ 2.0)) * (2.0 * x))"

// d/dx (x * exp(x)) = exp(x) + x*exp(x)   (product rule)
derive(Product(x, Exp(x)), x).toString
// res2: String = "(exp(x) + (x * exp(x)))"

Combine with simplify or simplifyFully to reduce the result:

simplifyFully(derive(Product(x, Exp(x)), x)).toString
// res3: String = "(exp(x) + (x * exp(x)))"

The hyperbolic and reciprocal-trigonometric functions (sinh/cosh/tanh, asinh/acosh/atanh, sec/csc/cot, sech/csch/coth) are first-class nodes with the standard rules:

// d/dx tanh(x) = sech(x)^2
derive(Tanh(x), x).toString
// res4: String = "(sech(x) ^ 2.0)"

Higher-order derivatives

Wrap the result in a second derive call (or use _Derivative nodes):

// d²/dx² x^4 = 12x²
val d2 = derive(derive(Power(x, _Number(4.0)), x), x)
// d2: _Expression = Product(
//   a = _Number(d = 4.0),
//   b = Product(
//     a = _Number(d = 3.0),
//     b = Power(base = _Variable(variable = "x"), exp = _Number(d = 2.0))
//   )
// )
simplifyFully(d2).toString
// res5: String = "(12.0 * (x ^ 2.0))"

Evaluating a derivative numerically

val slope = derive(Sin(x), x)   // should be cos(x)
// slope: _Expression = Cos(e = _Variable(variable = "x"))
val envPi2 = new Environment(5, Map("x" -> _Number(math.Pi / 2)))
// envPi2: Environment = it.grypho.scala.leonardo.core.Environment@7a9e656f
slope.eval(envPi2)               // cos(π/2) ≈ 0
// res6: Either[_Expression, _Value] = Right(
//   value = _Number(d = 6.123233995736766E-17)
// )

Indefinite integration

integrate(e, v) applies the symbolic rule table (linearity, power rule, exp/sin/cos, 1/x → log, linear-argument chain rule):

integrate(Power(x, _Number(2.0)), x).toString
// res7: String = "((x ^ 3.0) / 3.0)"

integrate(Sin(x), x).toString
// res8: String = "((-1.0 * cos(x)) / 1.0)"

// ∫ 1/x dx = log(x)
integrate(Ratio(_Number(1.0), x), x).toString
// res9: String = "ln(x)"

// Chain rule: ∫ sin(3x) dx = -cos(3x)/3
integrate(Sin(Product(_Number(3.0), x)), x).toString
// res10: String = "((-1.0 * cos((3.0 * x))) / 3.0)"

Beyond the linear chain rule, non-linear u-substitution closes ∫ f(g(x))·g'(x) dx: the engine tries candidate inner functions g, divides the integrand by g', and integrates in g only when the quotient is free of x.

// ∫ x·e^(x²) dx = e^(x²)/2   (u = x²)
integrate(Product(x, Exp(Power(x, _Number(2.0)))), x).toString
// res11: String = "(exp((x ^ 2.0)) / 2.0)"

Radical integrands close by trigonometric / hyperbolic substitution; √(a²+x²) and √(x²−a²) use the hyperbolic substitution, so the result is written with the asinh/acosh functions:

// ∫ dx/√(x²+1) = asinh(x)
integrate(Ratio(_Number(1.0), Power(Sum(Power(x, _Number(2.0)), _Number(1.0)), _Number(0.5))), x).toString
// res12: String = "asinh(x)"

Classic non-elementary integrals are answered with their named special functions — the trigonometric integrals Si/Ci, the exponential integral Ei, the logarithmic integral li, the Fresnel integrals fresnelS/fresnelC, and the error function erf — symbolic nodes with numeric kernels, so the antiderivative still evaluates:

// ∫ sin(x)/x dx = Si(x)
integrate(Ratio(Sin(x), x), x).toString
// res13: String = "Si(x)"

Unsupported forms are left as _Integral nodes (symbolic, not an error):

integrate(Sin(Product(x, x)), x).toString    // sin(x²) has no closed form
// res14: String = "integral(sin((x * x)), x)"

Vector calculus

The gradient grad, divergence div, curl curl, Laplacian laplacian, Jacobian jacobian and Hessian hessian take a scalar or vector field followed by the ordered coordinate tuple — the order is never inferred, because it fixes the order of the result’s components. A vector field is an n×1 matrix.

import it.grypho.scala.leonardo.vector.*
import it.grypho.scala.leonardo.matrix._Matrix

val y = _Variable("y")
// y: _Variable = _Variable(variable = "y")
// grad(x^2 * y) = [2xy, x^2]^T
_Grad(Product(Power(x, _Number(2.0)), y), Vector(x, y)).eval(env).toExpression.toString
// res15: String = "[[((2.0 * x) * y)], [(x ^ 2.0)]]"

The coordinate-free identities hold, which makes them the natural property tests: curl(grad f) = 0, div(curl F) = 0, and laplacian is defined as div ∘ grad so the two cannot disagree. Shapes that have no meaning — a curl outside three dimensions, a component/coordinate count mismatch, a repeated coordinate — stay symbolic rather than being guessed.

Cartesian, cylindrical and spherical coordinates are all supported, through one set of orthogonal-curvilinear formulas parameterised by the system’s scale factors — Cartesian is simply the case where they are all 1:

val r  = _Variable("r")
val th = _Variable("t")
val ph = _Variable("p")
val atPoint = new Environment(variables =
  Map("r" -> _Number(2.0), "t" -> _Number(0.7), "p" -> _Number(0.4)))

The Newtonian potential is harmonic away from the origin — ∇²(1/r) = 0 — and that is a sharp check on the scale factors, since a single wrong one breaks it. The symbolic form does not visibly collapse (simplify does no common-factor cancellation), so evaluate it:

_Laplacian(Ratio(_Number(1.0), r), Vector(r, th, ph), CoordinateSystem.Spherical).eval(atPoint)
// res16: Either[_Expression, _Value] = Right(value = _Number(d = 0.0))

cylindrical is (r, θ, z); spherical is (r, θ, φ) with θ the polar angle (physics) and sphericalmaths has θ azimuthal and φ polar (mathematics). All are three-dimensional, and the coordinates are identified by position, not by name — so the two spherical keywords differ in argument order, not in naming. Exchanging the last two coordinates flips the handedness of the basis, so curl, a pseudo-vector, carries the matching sign; the two conventions agree on the curl of the same physical field.

Definite integration (Simpson’s rule)

_DefIntegral(e, v, lo, hi) computes the definite integral numerically using composite Simpson’s rule. It uses a compiled Double ⇒ Double closure when the integrand is free of unresolvable nodes — no per-step allocation:

// ∫₀¹ x² dx = 1/3
_DefIntegral(Power(x, _Number(2.0)), x, _Number(0.0), _Number(1.0)).eval(env)
// res17: Either[_Expression, _Value] = Right(
//   value = _Number(d = 0.33333333333333315)
// )

// ∫₀π sin(x) dx = 2
_DefIntegral(Sin(x), x, _Number(0.0), _Number(math.Pi)).eval(env)
// res18: Either[_Expression, _Value] = Right(
//   value = _Number(d = 2.0000000000010805)
// )

// ∫₁ᵉ 1/x dx = 1  (ln e − ln 1)
_DefIntegral(Ratio(_Number(1.0), x), x, _Number(1.0), _Number(math.E)).eval(env)
// res19: Either[_Expression, _Value] = Right(
//   value = _Number(d = 1.000000000000286)
// )

Function sampling

sample(e, v, lo, hi, n, env) evaluates an expression over a uniform grid of n points in [lo, hi], returning Vector[(Double, Double)] with non-finite results silently dropped:

val pts = sample(Sin(x), x, 0.0, math.Pi, 5, env)

pts.length
// res20: Int = 5

pts.map { (xi: Double, yi: Double) => f"($xi%.4f, $yi%.4f)" }
// res21: Vector[String] = Vector(
//   "(0.0000, 0.0000)",
//   "(0.7854, 0.7071)",
//   "(1.5708, 1.0000)",
//   "(2.3562, 0.7071)",
//   "(3.1416, 0.0000)"
// )

The fast path compiles the expression to a Double ⇒ Double closure (no per-step allocation). A fallback per-step eval handles non-compilable nodes like _Derivative.

The Syntax extension gives method-call form:

import it.grypho.scala.leonardo.scalar.Syntax.*
val pts2 = Sin(x).sample(x, 0.0, math.Pi, 5, env)

pts2 == pts
// res22: Boolean = true