diff --git a/CHANGELOG.md b/CHANGELOG.md index 6805e411..a6d7dbbb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,87 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- **`Differential` operator** (spec 068): canonical SciML/ModelingToolkit-style + partial-differentiation operator. `Differential(x)(f)` for a declared function + `f(x, y)` returns a `DerivativeExpr` representing `∂f/∂x`. Composition handles + cross and higher-order partials: `Differential(y)(Differential(x)(f))` is + `∂²f/∂y∂x` (right-to-left Leibniz convention — see *Fixed* below); + `Differential(x)(Differential(x)(f))` is `∂²f/∂x²` (adjacent same-variable + steps collapse). Mono-variable functions also work: `Differential(t)(u)`. +- **Symbolics-compatible algebraic surface on `Differential`**: + `Differential(t)^n` for `n ≥ 0` returns an order-`n` operator — + `Differential(t)^2 === Differential(t, 2)` and `Differential(t)^0 === identity`, + matching Symbolics.jl. `Differential * Differential` composes via Julia's + generic `∘`, so `(Differential(y) * Differential(x))(f) == + (Differential(y) ∘ Differential(x))(f) == Differential(y)(Differential(x)(f))`. +- **`expand_derivatives`** (Symbolics.jl parity): forces evaluation of a + `DerivativeExpr` to a plain `GiacExpr` (delegates to GIAC's `diff`), and is + the identity on any other value. Lets code written in the Symbolics.jl style + (`expand_derivatives(Differential(x)(expr))`) run unchanged on Giac.jl. Note + that Giac.jl's bare-expression branch already evaluates eagerly, so the + no-op fallback covers the common case. +- **Two-argument `Differential(var, n)` shorthand**: `Differential(x, 2)(f)` + is exactly equivalent to `Differential(x)(Differential(x)(f))`, matching + Symbolics.jl's `Differential(x, 2)` form. The default `Differential(var)` + is `Differential(var, 1)`. Applies uniformly to function-form and + bare-expression operands. +- **Bare-expression differentiation** via `Differential` (spec 068): + `Differential(x)(x^2 + y*x)` returns the plain `GiacExpr` `2*x + y`, + matching SymPy.jl's `diff(expr, var)` and Symbolics.jl's `Differential`. +- **Multi-variable partial derivative support**: `DerivativeExpr` now records + an ordered sequence of `(variable, order)` differentiation steps, enabling + cross partials and arbitrary-depth composition. The single-step case is the + original mono-variable behavior. +- **Math-convention `n`-th derivative display**: high-order derivatives + (`n ≥ 4`) now render with parenthesized superscript notation in `Base.show` + rather than long strings of primes. So `D(u, 5)` displays as `D: u⁽⁵⁾(t)` + (instead of `D: u'''''(t)`). The same applies to `DerivativePoint` display. + `Base.string` for `DerivativePoint` keeps prime notation regardless of order + because GIAC consumes prime-form initial-condition strings. +- See `docs/migration/d_to_differential.md` for a full mapping table from the + deprecated `D` shapes to `Differential`. +- **`examples/07_odes_pdes.jl`**: new Pluto notebook walking through symbolic + ODE solutions (first-order, harmonic, damped, third-order, RLC, forced) and + PDE expressions (heat, wave, transport, separation of variables for the heat + equation) using the canonical `Differential` operator. Linked from + `docs/src/pluto.md`. + +### Deprecated + +- **`D` operator** (spec 068): every call shape (`D(u)`, `D(u, n)`, `D(D(u))`, + `D(d::DerivativeExpr, n)`) now emits a one-time-per-call-site + `Base.depwarn` pointing to the `Differential` replacement. `D` will be + removed in the next published release. The multi-arg ambiguity case + (`D(f)` on `f(x, y)`) raises `ArgumentError` instead of warning, since + the previous silent-default-to-first-variable behavior was a correctness + hazard. The `D(d::DerivativeExpr)` chained call also now requires the + underlying derivative to be mono-variable; partial derivatives of mixed + variables must be composed via `Differential(var)(d)`. + +### Changed + +- **`DerivativeExpr` internal layout** (spec 068): the `varname::String` and + `order::Int` fields are replaced by `steps::Vector{Tuple{String, Int}}`. + This is an internal change — no user code is expected to construct + `DerivativeExpr` directly. Public arithmetic (`+`, `-`, `*`, `/`, `^`), + equation (`~`), and string-conversion behavior is preserved. + +### Fixed + +- **Partial-derivative pretty-print convention (right-to-left Leibniz)**: the + multi-step `Base.show` for `DerivativeExpr` previously printed + `Differential(y)(Differential(x)(f))` as `D: ∂²f/∂x∂y`, which is consistent + with the index/`f_{xy}` convention but inconsistent with the operator-product + reading `(∂/∂y)(∂/∂x)f = ∂²f/∂y∂x`. The display now uses the right-to-left + Leibniz convention — the **rightmost** variable in `∂ⁿf/∂v₁…∂vₙ` is applied + **first**, the leftmost last — so the same expression now prints as + `D: ∂²f/∂y∂x`. The `steps` field (innermost-first) and the GIAC-side + `diff(diff(f(x,y),x),y)` string are unchanged; only the human-facing + ∂-notation flips. By Schwarz/Clairaut the underlying value is symmetric + for sufficiently smooth functions, so this is purely a notation fix. + ## [0.14.1] - 2026-05-11 ### Added diff --git a/docs/make.jl b/docs/make.jl index d4807ced..c9599d25 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -50,6 +50,9 @@ makedocs( ], "Held Commands" => "held_commands.md", "Tables.jl Compatibility" => "tables.md", + "Migration" => [ + "D → Differential" => "migration/d_to_differential.md", + ], "Extensions" => [ "Symbolics.jl" => "extensions/symbolics.md", "MathJSON.jl" => "extensions/mathjson.md", diff --git a/docs/src/mathematics/calculus.md b/docs/src/mathematics/calculus.md index 0ad690a9..24df9def 100644 --- a/docs/src/mathematics/calculus.md +++ b/docs/src/mathematics/calculus.md @@ -69,6 +69,63 @@ diff(x^2 * y^3, y) # Output: 3*x^2*y^2 ``` +### Using the `Differential` operator + +Beyond the function-form `Giac.Commands.diff`, Giac.jl also exposes a SciML/ModelingToolkit-style operator: `Differential(var)` is a callable that captures a single differentiation variable and applies partial differentiation to its operand. This pairs naturally with declared functions (`@giac_var f(x, y)`) and with bare symbolic expressions. + +```julia +using Giac + +@giac_var x y t f(x, y) u(t) + +# Bare-expression form (matches SymPy.jl `diff(expr, var)`) +Differential(x)(x^2 + y*x) # Returns: 2*x + y +Differential(x)(sin(x) * exp(x)) # Returns: cos(x)*exp(x) + sin(x)*exp(x) + +# Function-form: declared function, returns a DerivativeExpr +Dx = Differential(x) +Dy = Differential(y) +Dx(f) # ∂f/∂x +Dy(Dx(f)) # ∂²f/∂y∂x (cross partial — see "Notation" below) +Dx(Dx(f)) # ∂²f/∂x² (adjacent same-var steps collapse) + +# Mono-variable functions work the same way +Dt = Differential(t) +Dt(u) # u'(t) +Dt(Dt(u)) # u''(t) + +# Two-argument shorthand for n-th order (matches Symbolics.jl's Differential(x, 2)) +Differential(x, 3)(f) # ∂³f/∂x³ — same as Dx(Dx(Dx(f))) +Differential(t, 2)(u) # u''(t) — same as Dt(Dt(u)) +Differential(x, 2)(x^3) # 6*x — bare-expression second derivative + +# Symbolics-style algebraic surface +Differential(t)^2 # === Differential(t, 2) +Differential(t)^0 # === identity +Differential(y) * Differential(x) # composition — same as Differential(y) ∘ Differential(x) +(Differential(y) ∘ Differential(x))(f) # ∂²f/∂y∂x + +# Force evaluation of a lazy DerivativeExpr (no-op on plain GiacExpr — +# matches Symbolics.jl's expand_derivatives) +expand_derivatives(Dx(f)) # → diff(f(x,y), x) as a GiacExpr +expand_derivatives(x^2 + y*x) # → x^2 + y*x (identity) +``` + +#### Notation: right-to-left Leibniz convention + +In partial-derivative output (`Base.show` and the math comments above), +Giac.jl uses the **right-to-left Leibniz convention**: in `∂ⁿf/∂v₁…∂vₙ`, +the **rightmost** variable in the denominator is applied **first**, the +**leftmost** last. This matches the operator product +`(∂/∂v₁)·…·(∂/∂vₙ) f` read as a chain of left-multiplications: +`Differential(y)(Differential(x)(f))` (apply `Differential(x)` first, then +`Differential(y)`) is written `∂²f/∂y∂x` and pretty-prints as +`D: ∂²f/∂y∂x`. By Schwarz/Clairaut the value is symmetric for +sufficiently smooth `f`, but the notation, the `Base.show` output, and the +internal `DerivativeExpr.steps` field track the order of application. + +The function-form path returns a `DerivativeExpr` that records the differentiation history as an ordered sequence of `(variable, order)` steps (innermost-first). The bare-expression path returns a plain `GiacExpr` (already simplified by GIAC). For high-order derivatives (`n ≥ 4`), the human-readable display switches from primes to math-convention parenthesized superscripts: `D(u, 5)` shows as `u⁽⁵⁾(t)` rather than `u'''''(t)`. See [Differential Equations](differential_equations.md) for ODE/PDE workflows and [the migration guide](../migration/d_to_differential.md) for porting from the deprecated `D` operator. + ## Integration ### Indefinite Integrals diff --git a/docs/src/mathematics/differential_equations.md b/docs/src/mathematics/differential_equations.md index b7382c74..87f55da4 100644 --- a/docs/src/mathematics/differential_equations.md +++ b/docs/src/mathematics/differential_equations.md @@ -1,44 +1,47 @@ # Differential Equations -Giac.jl provides symbolic solving of ordinary differential equations (ODEs) using GIAC's `desolve` command. The package includes a `D` operator following SciML/ModelingToolkit conventions for expressing derivatives naturally. +Giac.jl provides symbolic solving of ordinary differential equations (ODEs) using GIAC's `desolve` command. The package includes a `Differential` operator following SciML/ModelingToolkit conventions for expressing derivatives naturally, plus a deprecated `D` alias kept during the transition window. -## The D Operator +## The `Differential` Operator -The `D` operator provides a clean, Julian syntax for expressing derivatives: +`Differential(var)` is the canonical, SciML-style partial-differentiation operator. It captures a single variable and is applied as a callable. For a mono-variable function declared with `@giac_var t u(t)`, `Differential(t)(u)` represents `u'(t)`. Repeated application gives higher orders; mixed variables compose for partial derivatives. See [Calculus](calculus.md) for a multi-variable walkthrough. ```julia using Giac using Giac.Commands: desolve @giac_var t u(t) +Dt = Differential(t) -# Create derivative expressions -D(u) # First derivative u' -D(D(u)) # Second derivative u'' (chained) -D(u, 2) # Second derivative u'' (direct) -D(u, 3) # Third derivative u''' +Dt(u) # u'(t) +Dt(Dt(u)) # u''(t) (adjacent same-var steps collapse) ``` -### Comparison with Raw Syntax +### Comparison with raw GIAC syntax -| D Operator | Raw GIAC | Description | -|------------|----------|-------------| -| `D(u)` | `diff(u, t)` | First derivative | -| `D(D(u))` | `diff(diff(u, t), t)` | Second derivative (chained) | -| `D(u, 2)` | `diff(u, t, 2)` | Second derivative (direct) | -| `D(u)(0) ~ 1` | `"u'(0)=1"` | Initial condition for u'(0) | +| `Differential` form | Raw GIAC | Description | +|---|---|---| +| `Differential(t)(u)` | `diff(u(t), t)` | First derivative | +| `Differential(t)(Differential(t)(u))` | `diff(u(t), t, 2)` | Second derivative | +| `Differential(t)(u)(0) ~ 1` | `"u'(0)=1"` | Initial condition for `u'(0)` | ### Unicode identifiers -Function and variable names may use any Unicode letters that Julia and GIAC -accept as identifiers — Greek letters, mathematical italics, etc.: +Function and variable names may use any Unicode letters that Julia and GIAC accept: ```julia @giac_var 𝑧 ϕ(𝑧) -D(ϕ) # diff(ϕ(𝑧), 𝑧) -D(ϕ, 2) # diff(ϕ(𝑧), 𝑧, 2) +Dz = Differential(𝑧) +Dz(ϕ) # diff(ϕ(𝑧), 𝑧) +Dz(Dz(ϕ)) # diff(ϕ(𝑧), 𝑧, 2) ``` +### The deprecated `D` alias + +Earlier versions of Giac.jl exposed a `D` operator (e.g. `D(u)`, `D(u, 2)`, `D(D(u))`). `D` is still available during a one-release deprecation window, with each call site emitting a `Base.depwarn` pointing to the canonical replacement. `D` will be removed in the next published release. See [`docs/migration/d_to_differential.md`](../migration/d_to_differential.md) for the full mapping table. + +The one `D` call shape that *no longer works* even with a warning: `D(f)` for a multi-argument function like `@giac_var x y f(x, y)` now raises `ArgumentError` instead of silently differentiating against the first variable. Use `Differential(x)(f)` (or `Differential(y)(f)`) explicitly. + ## First-Order ODEs ### Basic Example @@ -50,9 +53,10 @@ using Giac using Giac.Commands: desolve @giac_var t u(t) tau U0 +Dt = Differential(t) # Define ODE: τu' + u = U₀ -ode = tau * D(u) + u ~ U0 +ode = tau * Dt(u) + u ~ U0 # Initial condition: u(0) = 1 initial = u(0) ~ 1 @@ -66,9 +70,10 @@ result = desolve([ode, initial], t, :u) ```julia @giac_var t V(t) R C Vs +Dt = Differential(t) # Capacitor voltage ODE: RC·V' + V = Vs -ode = R * C * D(V) + V ~ Vs +ode = R * C * Dt(V) + V ~ Vs initial = V(0) ~ 0 result = desolve([ode, initial], t, :V) @@ -86,37 +91,30 @@ using Giac using Giac.Commands: desolve @giac_var t u(t) +Dt = Differential(t) -# Define ODE using chained D -ode = D(D(u)) + u ~ 0 +# Define ODE +ode = Dt(Dt(u)) + u ~ 0 # Initial conditions -u0 = u(0) ~ 1 # u(0) = 1 -du0 = D(u)(0) ~ 0 # u'(0) = 0 +u0 = u(0) ~ 1 # u(0) = 1 +du0 = Dt(u)(0) ~ 0 # u'(0) = 0 # Solve result = desolve([ode, u0, du0], t, :u) # Returns: cos(t) ``` -### Alternative Syntax with D(u, 2) - -```julia -# Same ODE using direct order specification -ode = D(u, 2) + u ~ 0 -result = desolve([ode, u(0) ~ 1, D(u)(0) ~ 0], t, :u) -# Returns: cos(t) -``` - ### Damped Oscillator Solve `u'' + 2ζω₀u' + ω₀²u = 0`: ```julia @giac_var t u(t) zeta omega0 +Dt = Differential(t) -ode = D(u, 2) + 2*zeta*omega0*D(u) + omega0^2*u ~ 0 -result = desolve([ode, u(0) ~ 1, D(u)(0) ~ 0], t, :u) +ode = Dt(Dt(u)) + 2*zeta*omega0*Dt(u) + omega0^2*u ~ 0 +result = desolve([ode, u(0) ~ 1, Dt(u)(0) ~ 0], t, :u) ``` ## Third-Order ODEs @@ -128,34 +126,36 @@ using Giac using Giac.Commands: desolve @giac_var t y(t) +Dt = Differential(t) # Define ODE -ode = D(y, 3) - y ~ 0 +ode = Dt(Dt(Dt(y))) - y ~ 0 # Initial conditions y0 = y(0) ~ 1 -dy0 = D(y)(0) ~ 1 -d2y0 = D(y, 2)(0) ~ 1 +dy0 = Dt(y)(0) ~ 1 +d2y0 = Dt(Dt(y))(0) ~ 1 # Solve result = desolve([ode, y0, dy0, d2y0], t, :y) # Returns: exp(t) ``` -## Using D in ODE Expressions +## Using `Differential` in ODE Expressions -The `D` operator supports arithmetic operations, making it natural to build ODE expressions: +The `Differential` operator's result type supports arithmetic, making it natural to build ODE expressions: ```julia @giac_var t u(t) a b c +Dt = Differential(t) # Build complex ODE expressions -ode1 = D(D(u)) + a*D(u) + b*u ~ c -ode2 = D(u, 2) - 4*D(u) + 4*u ~ 0 +ode1 = Dt(Dt(u)) + a*Dt(u) + b*u ~ c +ode2 = Dt(Dt(u)) - 4*Dt(u) + 4*u ~ 0 # Combine with other GiacExpr forcing = sin(t) -ode3 = D(D(u)) + u ~ forcing +ode3 = Dt(Dt(u)) + u ~ forcing ``` ## Important Notes @@ -172,24 +172,28 @@ desolve([ode, u(0) ~ 1], t, :u) desolve([ode, u(0) ~ 1], t, u) # May not work as expected ``` -### Initial Conditions with D +### Initial Conditions with `Differential` -The `D(u)(0)` syntax creates an unevaluated derivative condition that GIAC interprets correctly: +The `Differential(t)(u)(0)` syntax creates an unevaluated derivative condition that GIAC interprets correctly: ```julia -D(u)(0) ~ 1 # Creates "u'(0)=1" for GIAC -D(u, 2)(0) ~ 0 # Creates "u''(0)=0" for GIAC +Dt = Differential(t) +Dt(u)(0) ~ 1 # Creates "u'(0)=1" for GIAC +Dt(Dt(u))(0) ~ 0 # Creates "u''(0)=0" for GIAC ``` +This works because the underlying `DerivativeExpr` is mono-variable. Point-evaluation syntax for partial derivatives of multi-variable functions (e.g. PDE boundary conditions) is not supported in this release; see spec 068 Open Questions for the planned path. + ### Systems of ODEs GIAC can solve systems of first-order ODEs: ```julia @giac_var t x(t) y(t) +Dt = Differential(t) # dx/dt = y, dy/dt = -x -sys = [D(x) ~ y, D(y) ~ -x] +sys = [Dt(x) ~ y, Dt(y) ~ -x] initial = [x(0) ~ 1, y(0) ~ 0] # Solve as a system (pass both variables) @@ -204,9 +208,10 @@ Model radioactive decay: `dN/dt = -λN` ```julia @giac_var t N(t) lambda N0 +Dt = Differential(t) # Decay equation -ode = D(N) + lambda * N ~ 0 +ode = Dt(N) + lambda * N ~ 0 initial = N(0) ~ N0 result = desolve([ode, initial], t, :N) @@ -219,8 +224,9 @@ Exponential growth model: `dP/dt = rP` ```julia @giac_var t P(t) r P0 +Dt = Differential(t) -ode = D(P) - r * P ~ 0 +ode = Dt(P) - r * P ~ 0 initial = P(0) ~ P0 result = desolve([ode, initial], t, :P) @@ -233,8 +239,9 @@ Temperature change: `dT/dt = -k(T - T_env)` ```julia @giac_var t T(t) k T_env T0 +Dt = Differential(t) -ode = D(T) + k * (T - T_env) ~ 0 +ode = Dt(T) + k * (T - T_env) ~ 0 initial = T(0) ~ T0 result = desolve([ode, initial], t, :T) @@ -249,8 +256,10 @@ result = desolve([ode, initial], t, :T) ## API Reference ```@docs +Differential D DerivativeExpr DerivativePoint DerivativeCondition +expand_derivatives ``` diff --git a/docs/src/migration/d_to_differential.md b/docs/src/migration/d_to_differential.md new file mode 100644 index 00000000..040de4e8 --- /dev/null +++ b/docs/src/migration/d_to_differential.md @@ -0,0 +1,138 @@ +# Migrating from `D` to `Differential` + +Spec 068 introduces `Differential` as the canonical, SciML/ModelingToolkit-style +partial-differentiation operator and deprecates the legacy `D` operator. This +guide maps every `D(...)` call shape to its `Differential(...)` replacement so +you can update existing notebooks, scripts, and downstream packages +mechanically. + +## Why migrate + +Three reasons: + +1. **Multi-variable support**. `D(f)` on a multi-argument function like + `@giac_var x y f(x, y)` used to silently differentiate with respect to the + first variable (`x`), producing a wrong answer that looked plausible. + Calling `D(f)` now raises `ArgumentError`. `Differential(x)(f)` makes the + choice explicit. +2. **Bare-expression operands**. `Differential(x)(x^2 + y*x)` returns + `2*x + y` directly. `D` only accepts function-form expressions like + `u(t)` — bare expressions raise an error. `Differential` matches what + SymPy.jl's `diff(expr, var)` and Symbolics.jl's `Differential` already do. +3. **Composition reads naturally**. `Differential(y)(Differential(x)(f))` + tracks the cross partial `∂²f/∂y∂x` (right-to-left Leibniz convention: + the rightmost variable is applied first — `Differential(x)` is applied + before `Differential(y)`); mixed orders + (`Differential(x)(Differential(y)(Differential(x)(f)))`) work the same way. + +## Timeline + +`D` is deprecated in `0.15.0` and will be removed in the next published +release after `0.15.0`. Each unique `D(...)` call site emits exactly one +`Base.depwarn` per Julia session pointing to the canonical replacement. + +## Side-by-side mapping + +For mono-variable functions (`@giac_var t u(t)`): + +| Deprecated `D` form | Canonical `Differential` form | +|---|---| +| `D(u)` | `Differential(t)(u)` | +| `D(u, 2)` | `Differential(t, 2)(u)` *or* `Differential(t)(Differential(t)(u))` | +| `D(D(u))` | `Differential(t, 2)(u)` *or* `Differential(t)(Differential(t)(u))` | +| `D(u, 3)` | `Differential(t, 3)(u)` *or* `Differential(t)(Differential(t)(Differential(t)(u)))` | +| `D(u)(0) ~ 1` | `Differential(t)(u)(0) ~ 1` (initial-condition syntax unchanged) | + +The two-argument shorthand `Differential(t, n)` matches Symbolics.jl's `Differential(x, 2)` form and is exactly equivalent to applying `Differential(t)` `n` times. Use whichever reads better in your code — they produce the same `DerivativeExpr` value. + +For multi-variable functions (`@giac_var x y f(x, y)`): + +| Was (silently broken or unsupported) | Now | +|---|---| +| `D(f)` (silently picked `x`!) | `Differential(x)(f)` (explicit) or `Differential(y)(f)` | +| Not expressible directly with `D` | `Differential(y)(Differential(x)(f))` for `∂²f/∂y∂x` | +| Not expressible directly with `D` | `Differential(x)(Differential(x)(f))` for `∂²f/∂x²` | +| `D(f, x)` (transitional alias, deprecated) | `Differential(x)(f)` | +| `D(f, x, 2)` (transitional alias, deprecated) | `Differential(x)(Differential(x)(f))` | + +For bare expressions: + +| Was (raised) | Now | +|---|---| +| `D(x^2 + y*x)` (raised "requires a function expression") | `Differential(x)(x^2 + y*x) → 2*x + y` | +| `D(sin(x) * exp(x))` (raised) | `Differential(x)(sin(x) * exp(x))` | + +## Initial conditions for ODEs + +The mono-variable initial-condition syntax is unchanged because `Differential` +returns the same `DerivativeExpr` value: + +```julia +using Giac +using Giac.Commands: desolve + +@giac_var t u(t) +Dt = Differential(t) + +ode = Dt(Dt(u)) + u ~ 0 +ic1 = u(0) ~ 1 +ic2 = Dt(u)(0) ~ 0 # u'(0) = 0 — produces "u'(0)=0" via DerivativePoint + +desolve([ode, ic1, ic2], t, :u) # → cos(t) +``` + +PDE-style point-evaluation (e.g., `Differential(x)(f)(0, 0) ~ 1`) is *not* +supported in this release and raises a clear error. See spec 068 Open +Questions for the planned path. + +## `Differential` collision with `Symbolics.jl` + +Symbolics.jl exports a type also called `Differential`. When both packages are +loaded, Julia warns about the ambiguity and requires you to qualify the name: + +```julia +using Giac +using Symbolics + +@giac_var x y f(x, y) +Giac.Differential(x)(f) # Giac's path +@variables a b +Symbolics.Differential(a) # Symbolics' path +``` + +`Giac.Differential` and `Symbolics.Differential` are independent types — there +is no implicit interop in this release. This is the standard Julia +ambiguity-resolution pattern, the same one used across the SciML stack +(`DifferentialEquations.solve` vs `DifferentialEquations.ODESolution.solve`, +etc.). A bridging extension can be added later if a concrete need surfaces. + +## Mechanical migration recipe + +1. **Search-and-replace** simple `D(u) → Differential(t)(u)` in your codebase + (the variable `t` is whatever single argument `u` was declared with). +2. For `D(u, 2)`, replace with `Differential(t)(Differential(t)(u))` or assign + `Dt = Differential(t)` once and write `Dt(Dt(u))`. +3. For `D(D(u))` chained calls, rewrite as `Dt(Dt(u))`. +4. For multi-variable functions where you previously called `D(f)`, decide + which variable you actually want and write `Differential(x)(f)` or + `Differential(y)(f)` explicitly. +5. Keep initial-condition syntax untouched: `D(u)(0) ~ 1` becomes + `Differential(t)(u)(0) ~ 1`. + +## Bare-expression migration (new capability) + +If you were manually shelling out to `Giac.Commands.diff(expr, var)` to +differentiate bare expressions, you can now write: + +```julia +using Giac +@giac_var x y + +# Before +result = Giac.Commands.diff(x^2 + y*x, x) + +# After +result = Differential(x)(x^2 + y*x) # → 2*x + y +``` + +Both still work; `Differential` is the recommended form going forward. diff --git a/docs/src/pluto.md b/docs/src/pluto.md index 4e123df8..8b62a85e 100644 --- a/docs/src/pluto.md +++ b/docs/src/pluto.md @@ -22,6 +22,24 @@ using Pluto Pluto.run(notebook="examples/02_latex.jl") ``` +## Symbolic ODEs and PDEs + +A worked-example notebook covering symbolic ODE solutions via `desolve` (first-order, harmonic and damped oscillators, third-order, RLC, forced) and partial-differential-equation expressions (heat, wave, transport, separation of variables for the heat equation) using the canonical `Differential` operator from spec 068. All cells are reactive — change a coefficient or initial condition and downstream cells update. + +```julia +using Pluto +Pluto.run(notebook="examples/07_odes_pdes.jl") +``` + +## Other example notebooks + +- `examples/01_basics.jl` — symbolic operations primer +- `examples/02_latex.jl` — LaTeX-rendering demo (referenced above) +- `examples/03_examples.jl` — comprehensive showcase by domain +- `examples/04_plotting.jl` — symbolic-to-plotting pipelines +- `examples/05_programming.jl` — multi-line GIAC programs +- `examples/06_symbolics_bridge.jl` — Giac.jl × Symbolics.jl interop + See screenshots: ![screencapture-pluto-notebook](assets/screencapture-pluto-notebook-latex_demo.png) diff --git a/examples/03_examples.jl b/examples/03_examples.jl index f2a3d208..b4d64ce5 100644 --- a/examples/03_examples.jl +++ b/examples/03_examples.jl @@ -1017,7 +1017,7 @@ version = "17.7.0+0" # ╔═╡ Cell order: # ╟─b0c1d2e3-f4a5-6789-bcde-200000000000 -# ╟─2715140c-b9b1-45ba-b0f9-7346e0473bb1 +# ╠═2715140c-b9b1-45ba-b0f9-7346e0473bb1 # ╠═b0c1d2e3-f4a5-6789-bcde-200000000001 # ╠═b0c1d2e3-f4a5-6789-bcde-200000000002 # ╟─b0c1d2e3-f4a5-6789-bcde-200000000010 diff --git a/examples/07_odes_pdes.jl b/examples/07_odes_pdes.jl new file mode 100644 index 00000000..4f1c16fb --- /dev/null +++ b/examples/07_odes_pdes.jl @@ -0,0 +1,383 @@ +### A Pluto.jl notebook ### +# v0.20.24 + +using Markdown +using InteractiveUtils + +# ╔═╡ affa8790-1fd6-4467-bfda-23be0b73a882 +begin + using Pkg + # Pkg.add(PackageSpec(url="https://github.com/s-celles/Giac.jl")) + Pkg.develop(PackageSpec(path="..")) +end + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000001 +begin + using Giac + using Giac.Commands: desolve, simplify +end + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000000 +md""" +# Symbolic ODEs and PDEs with Giac.jl + +Closed-form solutions of ordinary differential equations through GIAC's `desolve` command, expressed with the SciML/ModelingToolkit-style `Differential` operator (spec 068). The same operator builds partial-derivative expressions for partial differential equations, even though GIAC's `desolve` itself targets ODEs only — for full numerical PDE solutions, hand the symbolic problem off to `MethodOfLines.jl` or another PDE solver. + +This notebook is reactive: change a coefficient or an initial condition and every downstream cell updates. + +--- + +## 1. The `Differential` operator in two minutes + +`Differential(var)` is a callable that captures a single differentiation variable. Apply it to any `GiacExpr` to obtain a partial derivative. + +- **Function-form** operands (declared via `@giac_var f(x, y)`) → returns a `DerivativeExpr` that knows its differentiation history. Supports `~` to build equations and `(d)(0)` to build initial conditions. +- **Bare-expression** operands → returns a plain `GiacExpr` (already simplified by GIAC), matching SymPy's `diff(expr, var)`. + +The two-argument shorthand `Differential(var, n)` matches Symbolics.jl's `Differential(x, 2)` form. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000002 +@giac_var t x y u(t) y_t(t) Q(t) X(x) T(t) ψ(x, t) θ(x, y) τ U₀ ω₀ ζ Ω L R C V_s c α λ + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000003 +md""" +For the rest of this notebook we'll alias the most common operator: +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000004 +Dt = Differential(t) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000005 +md""" +First and second time-derivatives of `u(t)`: +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000006 +Dt(u) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000007 +Dt(Dt(u)) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000008 +md""" +Or equivalently with the n-th-order shorthand: +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000009 +Differential(t, 2)(u) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-70000000000a +md""" +--- + +## 2. First-order linear ODE + +**Problem.** Solve `τ u'(t) + u(t) = U₀` with `u(0) = 1`. This is the canonical first-order RC-charging equation. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000010 +ode_1 = τ * Dt(u) + u ~ U₀ + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000011 +ic_1 = u(0) ~ 1 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000012 +sol_1 = desolve([ode_1, ic_1], t, :u) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000013 +md""" +The closed form: `u(t) = U₀ + (1 - U₀) · exp(-t/τ)` — what you'd derive by hand. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000014 +md""" +--- + +## 3. Harmonic oscillator (second-order, undamped) + +**Problem.** `u''(t) + u(t) = 0` with `u(0) = 1` and `u'(0) = 0`. Initial conditions for derivatives use the same `Differential(t)(u)(0)` syntax — internally a `DerivativePoint`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000020 +ode_2 = Dt(Dt(u)) + ω₀^2*u ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000021 +ic_2_pos = u(0) ~ 1 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000022 +ic_2_vel = Dt(u)(0) ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000023 +sol_2 = desolve([ode_2, ic_2_pos, ic_2_vel], t, :u) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000024 +md""" +Result: `cos(t)`. Try changing `ic_2_vel` to `Dt(u)(0) ~ 1` and the solution becomes `sin(t) + cos(t)` (or the equivalent canonical form GIAC chooses). +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000025 +md""" +--- + +## 4. Damped oscillator (second-order, with friction) + +**Problem.** `u''(t) + 2 ζ ω₀ u'(t) + ω₀² u(t) = 0` — the standard damped-oscillator form. The qualitative behavior depends on `ζ` (under-damped, critically-damped, over-damped). We solve the symbolic case with arbitrary `ζ` and `ω₀`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000030 +ode_3 = Dt(Dt(u)) + 2*ζ*ω₀*Dt(u) + ω₀^2 * u ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000031 +sol_3 = desolve([ode_3, u(0) ~ 1, Dt(u)(0) ~ 0], t, :u) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000032 +md""" +GIAC returns the closed form in terms of `ζ` and `ω₀`. The result is large but exact — and you can simplify or substitute numeric values to see specific regimes. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000033 +md""" +--- + +## 5. Third-order ODE + +**Problem.** `y'''(t) - y(t) = 0` with `y(0) = 1`, `y'(0) = 1`, `y''(0) = 1`. The `Differential(t, n)` shorthand keeps the cell readable for higher orders. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000040 +ode_4 = Differential(t, 3)(y_t) - y_t ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000041 +sol_4 = desolve( + [ode_4, + y_t(0) ~ 1, + Dt(y_t)(0) ~ 1, + Dt(Dt(y_t))(0) ~ 1], + t, :y_t, +) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000042 +md""" +The unique solution is `exp(t)`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000043 +md""" +--- + +## 6. RLC circuit (driven, second-order) + +**Problem.** A series RLC circuit driven by a constant source: `L Q''(t) + R Q'(t) + Q(t)/C = V_s`. Solve for the charge `Q(t)` with `Q(0) = 0` and `I(0) = Q'(0) = 0`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000051 +ode_5 = L*Dt(Dt(Q)) + R*Dt(Q) + Q/C ~ V_s + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000052 +sol_5 = desolve([ode_5, Q(0) ~ 0, Dt(Q)(0) ~ 0], t, :Q) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000053 +md""" +The result is symbolic in `L`, `R`, `C`, `V_s`. Substitute numeric values for any parameter to inspect the under/over/critically-damped regime — exactly as in §4. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000054 +md""" +--- + +## 7. Forced oscillator (non-homogeneous second-order ODE) + +**Problem.** A harmonic oscillator with a sinusoidal driving force at angular frequency `Ω`: `u''(t) + ω₀² u(t) = sin(Ω t)`. The closed form has a particular solution that grows linearly in `t` exactly when `Ω = ω₀` (resonance) — `desolve` returns the symbolic answer in `Ω` and `ω₀` so you can compare regimes. +""" + +# ╔═╡ fc3e0b7a-bb9b-4e35-ad05-d3b929daa200 +ode_6 = [Dt(Dt(u)) + ω₀^2 * u ~ sin(Ω * t), u(0) ~ 0, Dt(u)(0) ~ 0] + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000062 +sol_forced = desolve( + ode_6, + t, :u, +) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000063 +md""" +The result is symbolic in `Ω` and `ω₀`. Substitute `Ω => ω₀` and the solution diverges — that's the resonance singularity. + +> **Note on systems of ODEs.** GIAC's `desolve` solves single ODEs (and systems through some alternative entry points outside this notebook). If you need a coupled-system solver with closed-form output across vector unknowns, [`Symbolics.jl`](https://github.com/JuliaSymbolics/Symbolics.jl) + [`MTK.jl`](https://github.com/SciML/ModelingToolkit.jl) is the better tool — and `to_symbolics` / `to_giac` let you round-trip parts of the problem through GIAC. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000064 +md""" +--- + +## 8. Building partial differential equations + +`Differential` composes for partial derivatives just as cleanly. Declare a function of multiple variables and chain operators: +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000071 +Dx = Differential(x) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000072 +md""" +**Heat equation** (1D): `∂ψ/∂t = α · ∂²ψ/∂x²`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000073 +heat_eq = Dt(ψ) ~ α * Differential(x, 2)(ψ) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000074 +md""" +**Wave equation** (1D): `∂²ψ/∂t² = c² · ∂²ψ/∂x²`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000075 +wave_eq = Differential(t, 2)(ψ) ~ c^2 * Differential(x, 2)(ψ) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000076 +md""" +**Steady-state 2D heat equation** (Laplace's equation for the temperature field `θ(x, y)` in a heat-conducting plate at thermodynamic equilibrium, no time dependence): `∂²θ/∂x² + ∂²θ/∂y² = 0`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000077 +laplace_eq = Differential(x, 2)(θ) + Differential(y, 2)(θ) ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000078 +md""" +**Mixed second-order partials.** `∂²ψ/∂x∂t` records two distinct steps; same-variable applications collapse: +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000079 +Dt(Dx(ψ)) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-70000000007a +Dx(Dx(ψ)) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-70000000007b +md""" +The left value is the cross partial `∂²ψ/∂x∂t` (two steps preserved). The right value is `∂²ψ/∂x²` (steps collapsed because both are on `x`). +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-70000000007c +md""" +--- + +## 9. What about solving PDEs symbolically? + +GIAC's `desolve` is an **ODE** solver — passing a PDE expression to it does *not* in general return a closed-form PDE solution. For closed-form PDE solutions, a few options: + +1. **Separation of variables by hand.** Substitute `ψ(x, t) = X(x) · T(t)`, equate the ratios, and solve the two resulting ODEs with `desolve`. The notebook below sketches this for the heat equation. +2. **Symbolics.jl + MethodOfLines.jl.** For numerical PDE solutions, convert the symbolic PDE expression to Symbolics.jl (`to_symbolics`) and discretize. +3. **Specific known forms.** GIAC's `pde_separate` (when available) can do separation of variables for selected PDEs. + +Below: the heat equation reduced to two ODEs by separation of variables. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000081 +md""" +Spatial part: `X''(x) + λ · X(x) = 0`. Pick `λ > 0` (we get sinusoidal modes). +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000082 +spatial_ode = Differential(x, 2)(X) + λ * X ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000083 +spatial_sol = desolve([spatial_ode, X(0) ~ 0, Dx(X)(0) ~ 1], x, :X) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000084 +md""" +Temporal part: `T'(t) + α · λ · T(t) = 0` — exponential decay with rate `α λ`. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000085 +temporal_ode = Dt(T) + α * λ * T ~ 0 + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000086 +temporal_sol = desolve([temporal_ode, T(0) ~ 1], t, :T) + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000087 +md""" +The full mode `ψ(x, t) = X(x) · T(t)` decays in time with rate `α λ` while oscillating in space — the standard heat-equation diffusion mode. A complete solution sums infinitely many such modes weighted by the boundary data. +""" + +# ╔═╡ c2d3e4f5-a6b7-8901-cdef-700000000090 +md""" +--- + +## 10. Migration note: `D` vs `Differential` + +The earlier `D` operator (e.g. `D(u)`, `D(u, 2)`, `D(D(u))`) still works during the deprecation window with a `Base.depwarn`. The mapping to `Differential` is mechanical: + +| Old `D` form | New `Differential` form | +|---|---| +| `D(u)` | `Differential(t)(u)` | +| `D(u, 2)` | `Differential(t, 2)(u)` *or* `Differential(t)(Differential(t)(u))` | +| `D(D(u))` | `Differential(t, 2)(u)` | +| `D(u)(0) ~ 1` | `Differential(t)(u)(0) ~ 1` | + +The full table — including the multi-variable cases that the old `D` could not handle — lives in the [migration guide](../docs/src/migration/d_to_differential.md). + +For a notebook that compares Giac.jl with Symbolics.jl side-by-side, see [`06_symbolics_bridge.jl`](06_symbolics_bridge.jl). +""" + +# ╔═╡ Cell order: +# ╟─affa8790-1fd6-4467-bfda-23be0b73a882 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000000 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000001 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000002 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000003 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000004 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000005 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000006 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000007 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000008 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000009 +# ╟─c2d3e4f5-a6b7-8901-cdef-70000000000a +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000010 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000011 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000012 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000013 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000014 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000020 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000021 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000022 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000023 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000024 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000025 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000030 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000031 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000032 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000033 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000040 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000041 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000042 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000043 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000051 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000052 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000053 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000054 +# ╠═fc3e0b7a-bb9b-4e35-ad05-d3b929daa200 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000062 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000063 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000064 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000071 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000072 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000073 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000074 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000075 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000076 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000077 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000078 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000079 +# ╠═c2d3e4f5-a6b7-8901-cdef-70000000007a +# ╟─c2d3e4f5-a6b7-8901-cdef-70000000007b +# ╟─c2d3e4f5-a6b7-8901-cdef-70000000007c +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000081 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000082 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000083 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000084 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000085 +# ╠═c2d3e4f5-a6b7-8901-cdef-700000000086 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000087 +# ╟─c2d3e4f5-a6b7-8901-cdef-700000000090 diff --git a/src/Giac.jl b/src/Giac.jl index 6fe18f0b..6276a5ef 100644 --- a/src/Giac.jl +++ b/src/Giac.jl @@ -77,6 +77,9 @@ include("tables.jl") include("substitute.jl") include("build_function.jl") +# Differential operator (068-multivar-d-operator) +include("differential.jl") + # GenTypes module - Scoped enum for GIAC types (041-scoped-type-enum) include("gen_types.jl") @@ -99,7 +102,12 @@ include("Constants.jl") export GiacExpr, GiacContext, GiacMatrix, GiacError, HelpResult, GiacInput # Derivative operator (035-derivative-operator) -export D, DerivativeExpr, DerivativePoint, DerivativeCondition +# Differential is the canonical SciML-style form (068-multivar-d-operator). +# D is retained as a deprecated alias and will be removed in the next published release. +export D, Differential, DerivativeExpr, DerivativePoint, DerivativeCondition + +# expand_derivatives — Symbolics.jl-compatible no-op / DerivativeExpr forcer. +export expand_derivatives # Core functions export giac_eval, to_julia, list_commands, help_count diff --git a/src/differential.jl b/src/differential.jl new file mode 100644 index 00000000..ae024f8e --- /dev/null +++ b/src/differential.jl @@ -0,0 +1,190 @@ +# Spec 068 — multivar D operator +# Differential operator (canonical, SciML-style) — applies partial differentiation +# to function-form expressions and to bare GiacExpr operands. + +""" + Differential(var::GiacExpr) + Differential(var::GiacExpr, order::Int) + +Canonical partial-differentiation operator, SciML/ModelingToolkit-style. + +`Differential(var)` captures a single differentiation variable; the resulting +callable applies a first-order partial derivative to its operand. The two-argument +form `Differential(var, n)` produces an operator that applies the `n`-th order +partial derivative in one step — equivalent to applying `Differential(var)` `n` +times — and matches Symbolics.jl's `Differential(x, n)` shape. + +# Two operand forms + +- **Function-form**: `Differential(x)(f)` for `@giac_var x y f(x,y)` returns a + `DerivativeExpr` representing `∂f/∂x`. Compose by chaining: + `Differential(y)(Differential(x)(f))` is `∂²f/∂y∂x` (see *Notation* + below). The order-2 shape `Differential(x, 2)(f)` produces the same result + as `Differential(x)(Differential(x)(f))` (steps collapse). +- **Bare-expression**: `Differential(x)(x^2 + y*x)` returns a plain `GiacExpr` + (already simplified by GIAC) — `2*x + y` in this case. The order-2 shape + `Differential(x, 2)(x^3)` returns `6*x` (after simplification). Aligns with + SymPy.jl's `diff(expr, var)` and Symbolics.jl's `Differential`. + +# Notation (right-to-left Leibniz) + +`Giac.jl` uses the right-to-left Leibniz convention for `∂ⁿf/∂v₁…∂vₙ`: +the **rightmost** variable is applied **first**, the **leftmost** last — +consistent with the operator product `(∂/∂v₁)·…·(∂/∂vₙ) f`. So +`Differential(y)(Differential(x)(f))` (apply `Differential(x)` first, then +`Differential(y)`) is written `∂²f/∂y∂x`, and pretty-printed as `D: ∂²f/∂y∂x`. +By Schwarz/Clairaut the value is symmetric for sufficiently smooth `f`, but +the notation, the `Base.show` output, and the `steps` field track the order +of application. + +# Composition + +Repeated application on the same variable collapses adjacent steps: +`Differential(x)(Differential(x)(f))` stores `[("x", 2)]` rather than +`[("x", 1), ("x", 1)]`. Mixed variables compose freely: +`Differential(y)(Differential(x)(f))` stores `[("x", 1), ("y", 1)]` +(innermost-first) and prints as `diff(diff(f(x,y),x),y)` for GIAC. + +The operators `^`, `*`, and `∘` are also supported for compact composition, +matching Symbolics.jl: `Differential(t)^2 === Differential(t, 2)`, +`Differential(y) * Differential(x) ≡ Differential(y) ∘ Differential(x)`, and +`Differential(t)^0 === identity`. + +# Examples + +```julia +using Giac +@giac_var x y t f(x, y) u(t) + +Differential(x)(f) # ∂f/∂x +Differential(y)(Differential(x)(f)) # ∂²f/∂y∂x (x first, y next) + +Differential(t)(u) # u'(t) (mono-variable case) + +# Order-n shorthand (Symbolics.jl-style) +Differential(t, 2)(u) # u''(t) — same as Differential(t)(Differential(t)(u)) +Differential(x, 3)(f) # ∂³f/∂x³ + +Differential(x)(x^2 + y*x) # 2*x + y (bare expression) +Differential(x, 2)(x^3) # 6*x (bare expression, second derivative) +``` + +# Name collision with Symbolics.jl + +`Symbolics.Differential` exists too. When both packages are loaded, qualify: +`Giac.Differential(x)(f)` vs `Symbolics.Differential(x)`. The two are +independent types — there is no implicit interop in this release. See +`docs/migration/d_to_differential.md` for the recipe. + +# See also +- [`D`](@ref): Deprecated mono-variable operator. `D` is being removed in the + next published release; migrate `D(u) → Differential(t)(u)`. +- [`DerivativeExpr`](@ref): The function-form result type. +""" +struct Differential + var::GiacExpr + order::Int + + function Differential(var::GiacExpr, order::Int = 1) + order >= 1 || throw(ArgumentError("Differential order must be ≥ 1, got: $order")) + return new(var, order) + end +end + +function Base.show(io::IO, D::Differential) + if D.order == 1 + print(io, "Differential(", D.var, ")") + else + print(io, "Differential(", D.var, ", ", D.order, ")") + end +end + +# Application on a GiacExpr: dispatches between function-form (returns +# DerivativeExpr) and bare-expression (returns plain GiacExpr via GIAC's diff). +function (D::Differential)(operand::GiacExpr) + parsed = _parse_function_args(string(operand)) + if parsed === nothing + # Bare-expression branch (spec 068 US4): delegate to GIAC's diff(expr, var, n). + # Returns a plain GiacExpr (already simplified by GIAC) — no DerivativeExpr + # wrapper because there is no function-name / multi-step accounting to track. + result = operand + for _ in 1:D.order + result = Commands.diff(result, D.var) + end + return result + end + funcname, _args = parsed + var_str = string(D.var) + return DerivativeExpr(operand, funcname, Tuple{String, Int}[(var_str, D.order)]) +end + +# Composition: append a step to an existing DerivativeExpr. +# `Differential(x, n)(d)` adds an n-order step against `x` to d.steps, +# collapsing if the previous step is on the same variable. +function (D::Differential)(d::DerivativeExpr) + var_str = string(D.var) + new_steps = _append_step(d.steps, (var_str, D.order)) + return DerivativeExpr(d.base_expr, d.funcname, new_steps) +end + +# Symbolics.jl-style algebraic surface on Differential operators. +# `Differential(var)^n` yields the n-th order operator (n == 0 → identity, as +# in Symbolics). Composition via `*` and `∘` (Julia's generic ∘ already does +# the right thing on callables, but we add an explicit method so `*` can +# delegate uniformly and dispatch is clear in stack traces). + +""" + Base.:^(D::Differential, n::Integer) -> Differential | typeof(identity) + +Return a `Differential` whose order is multiplied by `n`. Mirrors +Symbolics.jl: `Differential(x)^n === Differential(x, D.order * n)` for `n ≥ 1`, +and `Differential(x)^0 === identity` (the identity function). + +`n` must be non-negative. Negative powers are not defined (anti-derivative +would require an explicit constant of integration). +""" +function Base.:^(D::Differential, n::Integer) + n >= 0 || throw(ArgumentError("Differential^n requires n ≥ 0, got: $n")) + n == 0 && return identity + return Differential(D.var, D.order * n) +end + +""" + Base.:*(D1::Differential, D2::Differential) -> ComposedFunction + +Compose two `Differential` operators: `(D1 * D2)(f) == D1(D2(f))`. Mirrors +Symbolics.jl's `Differential * Differential` and falls back to Julia's +generic `ComposedFunction` so applying the result on a `GiacExpr` or +`DerivativeExpr` reuses the existing `Differential` call-site logic +(including same-variable step collapse). +""" +Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2 + +""" + expand_derivatives(x) -> x + +Force evaluation of any lazy derivative wrapper, returning a plain `GiacExpr`. + +For `Differential(var)(f)` where `f` is a declared user function (`DerivativeExpr`), +`expand_derivatives` calls GIAC's `diff` and returns the resulting `GiacExpr`. For +everything else (already a `GiacExpr` produced by Giac.jl's eager bare-expression +branch, or any value not wrapped in `DerivativeExpr`), it is the identity — so +code written in the Symbolics.jl style (`expand_derivatives(Differential(x)(expr))`) +stays correct on Giac.jl without behavioral change. + +# Examples + +```julia +@giac_var x y f(x, y) + +expand_derivatives(x^2 + y) # identity → x^2+y +expand_derivatives(Differential(x)(x^2)) # identity → 2*x (already a GiacExpr) +expand_derivatives(Differential(x)(f)) # → diff(f(x,y),x) as a GiacExpr +``` + +# See also +- [`Differential`](@ref): The operator producing lazy or eager derivatives. +- [`DerivativeExpr`](@ref): The lazy function-form derivative wrapper. +""" +expand_derivatives(x) = x +expand_derivatives(d::DerivativeExpr) = _to_giac_expr(d) diff --git a/src/types.jl b/src/types.jl index 0f3fb946..f2754395 100644 --- a/src/types.jl +++ b/src/types.jl @@ -333,42 +333,87 @@ GiacExpr: y+2*x # Derivative Operator D (035-derivative-operator) # ============================================================================ +# Names treated as GIAC operations rather than user-declared functions. +# Shared between _parse_function_args and _parse_function_expr. +const _GIAC_OPERATIONS = Set([ + "diff", "integrate", "limit", "sum", "product", + "solve", "desolve", "simplify", "factor", "expand", + "sin", "cos", "tan", "exp", "log", "sqrt", "abs", +]) + """ - _parse_function_expr(expr_str::String) -> Union{Tuple{String, String}, Nothing} + _parse_function_args(expr_str::String) -> Union{Tuple{String, Vector{String}}, Nothing} -Parse a function expression to extract the function name and first variable. +Parse a function expression to extract the function name and the full list +of arguments. Strict generalization of [`_parse_function_expr`](@ref). -For expressions like "u(t)" returns ("u", "t"). -For expressions like "f(x,y)" returns ("f", "x"). -For non-function expressions, returns nothing. +Returns `nothing` for non-function expressions (bare symbols, composite +expressions, GIAC built-ins from `_GIAC_OPERATIONS`). # Examples ```julia -_parse_function_expr("u(t)") # ("u", "t") -_parse_function_expr("f(x,y)") # ("f", "x") -_parse_function_expr("x") # nothing +_parse_function_args("u(t)") # ("u", ["t"]) +_parse_function_args("f(x,y)") # ("f", ["x", "y"]) +_parse_function_args("g(a, b, c)") # ("g", ["a", "b", "c"]) +_parse_function_args("x") # nothing +_parse_function_args("sin(x)") # nothing (GIAC operation) ``` """ -function _parse_function_expr(expr_str::String) +function _parse_function_args(expr_str::String) # \p{L}/\p{N}: GIAC accepts Unicode identifiers (ϕ, 𝑧, α, …); # ASCII-only [a-zA-Z] would reject them. m = match(r"^([\p{L}_][\p{L}\p{N}_]*)\(([^)]+)\)$", expr_str) if m !== nothing funcname = m.captures[1] args_str = m.captures[2] - # Get first variable (strip whitespace) - varname = strip(split(args_str, ",")[1]) - # Exclude GIAC operations - giac_operations = Set(["diff", "integrate", "limit", "sum", "product", - "solve", "desolve", "simplify", "factor", "expand", - "sin", "cos", "tan", "exp", "log", "sqrt", "abs"]) - if funcname ∉ giac_operations - return (funcname, String(varname)) + if funcname ∉ _GIAC_OPERATIONS + args = String.(strip.(split(args_str, ","))) + return (funcname, args) end end return nothing end +""" + _parse_function_expr(expr_str::String) -> Union{Tuple{String, String}, Nothing} + +Parse a function expression to extract the function name and first variable. +Thin wrapper around [`_parse_function_args`](@ref) that projects to the first +argument — preserved for backward compatibility with call sites that only +need the first variable. + +For expressions like "u(t)" returns ("u", "t"). +For expressions like "f(x,y)" returns ("f", "x"). +For non-function expressions, returns nothing. +""" +function _parse_function_expr(expr_str::String) + parsed = _parse_function_args(expr_str) + parsed === nothing && return nothing + funcname, args = parsed + return (funcname, args[1]) +end + +""" + _append_step(steps, step) -> Vector{Tuple{String, Int}} + +Append a `(var, order)` differentiation step to a `Vector{Tuple{String, Int}}`, +collapsing the new step into the previous one when they share a variable. + +Used by both `D(::DerivativeExpr, ...)` and `Differential` composition to keep +multi-step differentiation history minimal: `[(x, 1)]` + `(x, 1)` → `[(x, 2)]`. +Non-adjacent same-variable steps are deliberately *not* collapsed — Schwarz/ +Clairaut commutativity is a math property, but the steps vector records what +the user wrote, in order. +""" +function _append_step(steps::Vector{Tuple{String, Int}}, step::Tuple{String, Int}) + if !isempty(steps) && last(steps)[1] == step[1] + last_var, last_n = last(steps) + return Tuple{String, Int}[steps[1:end-1]..., (last_var, last_n + step[2])] + else + return Tuple{String, Int}[steps..., step] + end +end + """ DerivativeExpr @@ -409,8 +454,7 @@ ode = D(D(u)) + u ~ 0 # Equivalent to diff(u,t,2) + u = 0 struct DerivativeExpr base_expr::GiacExpr funcname::String - varname::String - order::Int + steps::Vector{Tuple{String, Int}} end """ @@ -461,49 +505,188 @@ desolve([ode, y(0) ~ 1, D(y)(0) ~ 1, D(y,2)(0) ~ 1], t, y) - `desolve`: Solving differential equations (via `Giac.Commands`) """ function D(expr::GiacExpr) - parsed = _parse_function_expr(string(expr)) + parsed = _parse_function_args(string(expr)) if parsed === nothing throw(ArgumentError("D() requires a function expression like u(t), got: $(string(expr))")) end - funcname, varname = parsed - return DerivativeExpr(expr, funcname, varname, 1) + funcname, args = parsed + if length(args) >= 2 + # Multi-arg ambiguity guard — raises, does NOT depwarn (the error itself + # carries the deprecation note). + throw(ArgumentError( + "$funcname has multiple arguments (" * join(args, ", ") * + "); specify which variable to differentiate with respect to. " * + "Use the canonical form: Differential(" * args[1] * ")(" * funcname * + ") or Differential(" * args[2] * ")(" * funcname * + "). Note: D is deprecated and will be removed in the next published release." + )) + end + Base.depwarn( + "D($funcname) is deprecated; use Differential($(args[1]))($funcname) instead. " * + "D will be removed in the next published release.", + :D, + ) + return DerivativeExpr(expr, funcname, Tuple{String, Int}[(args[1], 1)]) end function D(expr::GiacExpr, n::Int) if n < 1 throw(ArgumentError("Derivative order must be positive, got: $n")) end - parsed = _parse_function_expr(string(expr)) + parsed = _parse_function_args(string(expr)) if parsed === nothing throw(ArgumentError("D() requires a function expression like u(t), got: $(string(expr))")) end - funcname, varname = parsed - return DerivativeExpr(expr, funcname, varname, n) + funcname, args = parsed + if length(args) >= 2 + throw(ArgumentError( + "$funcname has multiple arguments (" * join(args, ", ") * + "); specify which variable to differentiate with respect to. " * + "Use the canonical form: Differential(" * args[1] * ")(" * funcname * + ") or Differential(" * args[2] * ")(" * funcname * + "). Note: D is deprecated and will be removed in the next published release." + )) + end + Base.depwarn( + "D($funcname, $n) is deprecated; use Differential($(args[1])) applied $n times to $funcname instead. " * + "D will be removed in the next published release.", + :D, + ) + return DerivativeExpr(expr, funcname, Tuple{String, Int}[(args[1], n)]) end function D(d::DerivativeExpr) - return DerivativeExpr(d.base_expr, d.funcname, d.varname, d.order + 1) + if length(d.steps) != 1 + throw(ArgumentError( + "D(::DerivativeExpr) without a variable is only supported for mono-variable derivatives. " * + "Use Differential(var)(d) instead." + )) + end + Base.depwarn( + "D(::DerivativeExpr) is deprecated; use Differential(var)(d) instead. " * + "D will be removed in the next published release.", + :D, + ) + var, n = d.steps[1] + return DerivativeExpr(d.base_expr, d.funcname, Tuple{String, Int}[(var, n + 1)]) end function D(d::DerivativeExpr, n::Int) if n < 1 throw(ArgumentError("Derivative order must be positive, got: $n")) end - return DerivativeExpr(d.base_expr, d.funcname, d.varname, d.order + n) + if length(d.steps) != 1 + throw(ArgumentError( + "D(::DerivativeExpr, ::Int) without a variable is only supported for mono-variable derivatives. " * + "Use Differential(var)(d) instead." + )) + end + Base.depwarn( + "D(::DerivativeExpr, $n) is deprecated; use Differential(var) applied $n times to d instead. " * + "D will be removed in the next published release.", + :D, + ) + var, k = d.steps[1] + return DerivativeExpr(d.base_expr, d.funcname, Tuple{String, Int}[(var, k + n)]) end -# String conversion - produces diff notation for use in equations -function Base.string(d::DerivativeExpr) - if d.order == 1 - return "diff($(string(d.base_expr)),$(d.varname))" +# spec 068 US5 — transitional D(f, x) alias overloads +# Each emits a depwarn pointing to Differential and delegates. Note: the +# Differential type is defined later in src/differential.jl, but these +# function bodies resolve `Differential` at call time, so forward reference +# is fine. + +function D(expr::GiacExpr, var::GiacExpr) + Base.depwarn( + "D(expr, var) is deprecated; use Differential($(string(var)))(expr) instead. " * + "D will be removed in the next published release.", + :D, + ) + return Differential(var)(expr) +end + +function D(expr::GiacExpr, var::GiacExpr, n::Int) + if n < 1 + throw(ArgumentError("Derivative order must be positive, got: $n")) + end + Base.depwarn( + "D(expr, var, $n) is deprecated; use Differential($(string(var))) applied $n times to expr instead. " * + "D will be removed in the next published release.", + :D, + ) + parsed = _parse_function_args(string(expr)) + if parsed === nothing + throw(ArgumentError("D(expr, var, n) requires a function expression like u(t), got: $(string(expr))")) + end + funcname, _args = parsed + return DerivativeExpr(expr, funcname, Tuple{String, Int}[(string(var), n)]) +end + +function D(d::DerivativeExpr, var::GiacExpr) + Base.depwarn( + "D(::DerivativeExpr, var) is deprecated; use Differential($(string(var)))(d) instead. " * + "D will be removed in the next published release.", + :D, + ) + return Differential(var)(d) +end + +# Helper for ∂-style display: superscript a non-negative integer (e.g. 12 → "¹²"). +function _superscript_int(n::Int) + digits_map = Dict('0'=>'⁰','1'=>'¹','2'=>'²','3'=>'³','4'=>'⁴','5'=>'⁵', + '6'=>'⁶','7'=>'⁷','8'=>'⁸','9'=>'⁹') + return join(digits_map[c] for c in string(n)) +end + +# Math-convention threshold: primes for n ≤ 3, parenthesized superscript ⁽ⁿ⁾ for n ≥ 4. +# Spec 068 — keeps the familiar u', u'', u''' for everyday ODE work and switches +# to u⁽⁴⁾(t), u⁽⁵⁾(t), … for higher orders. +const _PRIMES_MAX = 3 + +# Render a derivative-order suffix on a function name in "show" output. +# n=1 → "'", n=2 → "''", n=3 → "'''", n=4 → "⁽⁴⁾", n=12 → "⁽¹²⁾", etc. +function _derivative_order_suffix(n::Int) + if n <= _PRIMES_MAX + return repeat("'", n) else - return "diff($(string(d.base_expr)),$(d.varname),$(d.order))" + return "⁽" * _superscript_int(n) * "⁾" end end +# String conversion - produces diff notation for use in equations. +# Mono-step or all-same-variable: emit a single diff(expr, var, n) call. +# Otherwise: emit nested diff(...) calls in the order steps were applied. +function Base.string(d::DerivativeExpr) + base_str = string(d.base_expr) + if length(d.steps) == 1 + var, n = d.steps[1] + return n == 1 ? "diff($base_str,$var)" : "diff($base_str,$var,$n)" + end + s = base_str + for (var, n) in d.steps + s = n == 1 ? "diff($s,$var)" : "diff($s,$var,$n)" + end + return s +end + function Base.show(io::IO, d::DerivativeExpr) - primes = repeat("'", d.order) - print(io, "D: ", d.funcname, primes, "(", d.varname, ")") + if length(d.steps) == 1 + var, n = d.steps[1] + suffix = _derivative_order_suffix(n) + print(io, "D: ", d.funcname, suffix, "(", var, ")") + else + total = sum(n for (_, n) in d.steps) + ord_sup = total == 1 ? "" : _superscript_int(total) + # Right-to-left Leibniz convention: ∂ⁿf/∂v₁…∂vₙ reads with v₁ + # last-applied (leftmost) and vₙ first-applied (rightmost). `d.steps` + # records innermost-first (application order), so reverse it. + denom = IOBuffer() + for (var, n) in Iterators.reverse(d.steps) + print(denom, "∂", var) + n == 1 || print(denom, _superscript_int(n)) + end + print(io, "D: ∂", ord_sup, d.funcname, "/", String(take!(denom))) + end end # ============================================================================ @@ -538,12 +721,17 @@ struct DerivativePoint end function Base.string(dp::DerivativePoint) + # GIAC parses prime notation for derivative initial conditions (e.g. u'(0)=1). + # We keep `'`-style here regardless of order so the string is GIAC-consumable; + # the human-readable display falls through to `Base.show` below, which uses + # the math-convention u⁽ⁿ⁾ form for n ≥ 4. primes = repeat("'", dp.order) return dp.funcname * primes * "(" * join(dp.point_args, ",") * ")" end function Base.show(io::IO, dp::DerivativePoint) - print(io, string(dp)) + suffix = _derivative_order_suffix(dp.order) + print(io, dp.funcname, suffix, "(", join(dp.point_args, ","), ")") end # Callable - produces DerivativePoint for initial conditions @@ -564,8 +752,15 @@ D(u, 2)(0) ~ 0 # Creates GiacExpr: "u''(0)=0" ``` """ function (d::DerivativeExpr)(args...) + if length(d.steps) >= 2 + throw(ArgumentError( + "Point evaluation is not supported for partial derivatives in this release. " * + "See spec 068 Open Questions for the planned PDE-initial-condition path." + )) + end + _, n = d.steps[1] arg_strs = [_arg_to_giac_string(arg) for arg in args] - return DerivativePoint(d.funcname, d.order, arg_strs) + return DerivativePoint(d.funcname, n, arg_strs) end """ diff --git a/test/runtests.jl b/test/runtests.jl index 3413f3c1..1f061f1c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,9 @@ using LinearAlgebra # Macro tests (011-giac-symbol-macro) include("test_macros.jl") + # Multivar derivative operator tests (068-multivar-d-operator) + include("test_differential.jl") + # Matrix display tests (011-giacmatrix-display) include("test_matrix_display.jl") diff --git a/test/test_differential.jl b/test/test_differential.jl new file mode 100644 index 00000000..e2fc4f30 --- /dev/null +++ b/test/test_differential.jl @@ -0,0 +1,346 @@ +# Spec 068 — multivar D operator +# Tests for Differential and the multi-variable derivative surface. +using Test +using Giac + +@testset "spec 068 — multivar D operator" begin + + # ------------------------------------------------------------------ + # Foundational helpers: _append_step + # ------------------------------------------------------------------ + @testset "_append_step collapses adjacent same-variable steps" begin + T = Tuple{String, Int} + + # empty + step → singleton + @test Giac._append_step(T[], ("x", 1)) == T[("x", 1)] + + # same-var adjacent → collapse + @test Giac._append_step(T[("x", 1)], ("x", 1)) == T[("x", 2)] + @test Giac._append_step(T[("x", 2)], ("x", 3)) == T[("x", 5)] + + # different var → append (no collapse) + @test Giac._append_step(T[("x", 1)], ("y", 1)) == T[("x", 1), ("y", 1)] + + # non-adjacent same-var → append (no collapse) + @test Giac._append_step(T[("x", 1), ("y", 1)], ("x", 1)) == + T[("x", 1), ("y", 1), ("x", 1)] + + # only the last step is collapsed against + @test Giac._append_step(T[("y", 1), ("x", 1)], ("x", 1)) == + T[("y", 1), ("x", 2)] + end + + # ------------------------------------------------------------------ + # US1 — canonical Differential, function-form mono-variable + # ------------------------------------------------------------------ + @testset "Differential mono-variable function-form (US1)" begin + @giac_var t u(t) + Dt = Differential(t) + du = Dt(u) + @test du isa DerivativeExpr + @test du.steps == [("t", 1)] + @test du.funcname == "u" + @test string(du) == "diff(u(t),t)" + + # Show on Differential prints variable + io = IOBuffer() + show(io, Dt) + @test String(take!(io)) == "Differential(t)" + end + + # ------------------------------------------------------------------ + # US1 — multi-variable function-form (partials, cross, higher-order) + # ------------------------------------------------------------------ + @testset "Differential multi-variable function-form (US1)" begin + @giac_var x y f(x, y) + Dx = Differential(x) + Dy = Differential(y) + + # Simple partials + df_dx = Dx(f) + @test df_dx isa DerivativeExpr + @test df_dx.steps == [("x", 1)] + @test df_dx.funcname == "f" + @test string(df_dx) == "diff(f(x,y),x)" + + df_dy = Dy(f) + @test df_dy.steps == [("y", 1)] + @test string(df_dy) == "diff(f(x,y),y)" + + # Higher-order (same var) — adjacent steps collapse + d2f_dx2 = Dx(Dx(f)) + @test d2f_dx2.steps == [("x", 2)] + @test string(d2f_dx2) == "diff(f(x,y),x,2)" + + # Cross partial — distinct vars don't collapse + d2f_dxdy = Dy(Dx(f)) + @test d2f_dxdy.steps == [("x", 1), ("y", 1)] + @test string(d2f_dxdy) == "diff(diff(f(x,y),x),y)" + end + + # ------------------------------------------------------------------ + # US1 — three-variable composition + value equivalence + # ------------------------------------------------------------------ + @testset "Differential three-variable composition + value (US1)" begin + @giac_var x y z g(x, y, z) + d_xyz = Differential(z)(Differential(y)(Differential(x)(g))) + @test d_xyz.steps == [("x", 1), ("y", 1), ("z", 1)] + @test string(d_xyz) == "diff(diff(diff(g(x,y,z),x),y),z)" + + # Value equivalence: ∂/∂x of (x^2*y + sin(y)) → 2*x*y + @giac_var a b + Dx = Differential(a) + # bare-expression form lands in US4 — for now, build through the function-form + # by declaring the expression as a function. + @giac_var h(a, b) + # We can't redefine h to equal a specific expression, but we can verify + # that string(Differential(a)(h)) is "diff(h(a,b),a)" — value-side check + # of the bare-expression flavor lands in US4. + @test string(Dx(h)) == "diff(h(a,b),a)" + end + + # ------------------------------------------------------------------ + # US4 — bare-expression differentiation (Differential matches SymPy/Symbolics) + # ------------------------------------------------------------------ + @testset "Differential on bare expressions (US4)" begin + @giac_var x y z + + # ∂/∂x of x^2 + y*x = 2x + y + result = Differential(x)(x^2 + y*x) + @test result isa GiacExpr + @test !(result isa DerivativeExpr) + # Compare via simplification: 2*x + y + s = string(result) + # GIAC may return "2*x+y" or "y+2*x" depending on canonicalization; check both terms. + @test occursin("2*x", s) && occursin("y", s) + + # ∂/∂x of sin(x) is cos(x) (round-trip through GIAC's diff) + result_sin = Differential(x)(sin(x)) + @test result_sin isa GiacExpr + @test occursin("cos", string(result_sin)) + + # ∂²/∂x² of x^3 = 6*x (GIAC may print as "3*2*x" without auto-simplifying; + # use simplify to canonicalize before comparing) + result_x3 = Differential(x)(Differential(x)(x^3)) + @test result_x3 isa GiacExpr + simplified_x3 = Giac.Commands.simplify(result_x3) + @test occursin("6*x", string(simplified_x3)) + + # Cross partial: ∂²/∂x∂y of x*y^2 = 2*y + result_xy = Differential(y)(Differential(x)(x*y^2)) + @test result_xy isa GiacExpr + simplified_xy = Giac.Commands.simplify(result_xy) + @test occursin("2*y", string(simplified_xy)) + + # Absent variable → 0 (CAS convention) + result_zero = Differential(z)(x^2 + y) + @test result_zero isa GiacExpr + @test string(result_zero) == "0" + end + + # ------------------------------------------------------------------ + # Differential(var, n) — Symbolics.jl-style n-th order constructor + # ------------------------------------------------------------------ + @testset "Differential(var, n) constructor (068-multivar-d-operator)" begin + @giac_var x y t f(x, y) u(t) + + @testset "constructor + show" begin + D1 = Differential(x) + @test D1.order == 1 + io = IOBuffer(); show(io, D1) + @test String(take!(io)) == "Differential(x)" + + D2 = Differential(x, 2) + @test D2.order == 2 + io = IOBuffer(); show(io, D2) + @test String(take!(io)) == "Differential(x, 2)" + + # invalid orders + @test_throws ArgumentError Differential(x, 0) + @test_throws ArgumentError Differential(x, -1) + end + + @testset "function-form: Differential(x, n)(f) ≡ Dx applied n times" begin + # Order 2 same-var + d_via_n = Differential(x, 2)(f) + d_via_chain = Differential(x)(Differential(x)(f)) + @test d_via_n.steps == d_via_chain.steps == [("x", 2)] + + # Order 3 mono-variable + d3 = Differential(t, 3)(u) + @test d3.steps == [("t", 3)] + @test string(d3) == "diff(u(t),t,3)" + end + + @testset "function-form: composition with Differential(_, n)" begin + # Differential(x, 2) on top of an existing Differential(y)(f) + mid = Differential(y)(f) # steps = [("y", 1)] + full = Differential(x, 2)(mid) # append ("x", 2) — distinct var, no collapse + @test full.steps == [("y", 1), ("x", 2)] + + # Same-variable: Differential(y, 2) on top of Differential(y)(f) collapses + mid2 = Differential(y)(f) # steps = [("y", 1)] + full2 = Differential(y, 2)(mid2) # append ("y", 2) — collapses to ("y", 3) + @test full2.steps == [("y", 3)] + end + + @testset "bare-expression: Differential(x, n)(expr)" begin + # Differential(x, 2)(x^3) = 6*x (after simplify) + result = Differential(x, 2)(x^3) + @test result isa GiacExpr + @test !(result isa DerivativeExpr) + simplified = Giac.Commands.simplify(result) + @test occursin("6*x", string(simplified)) + end + end + + @testset "Return-type discipline: function-form vs bare (US4)" begin + @giac_var x y f(x, y) + # function-form → DerivativeExpr + @test Differential(x)(f) isa DerivativeExpr + # bare-expression → plain GiacExpr (NOT DerivativeExpr) + @test Differential(x)(x^2 + y*x) isa GiacExpr + @test !(Differential(x)(x^2 + y*x) isa DerivativeExpr) + end + + # ------------------------------------------------------------------ + # US1 — ODE workflow (Differential composes with desolve) + # ------------------------------------------------------------------ + @testset "Differential ODE workflow with desolve (US1)" begin + using Giac.Commands: desolve + @giac_var t u(t) + Dt = Differential(t) + + # u'' + u = 0, u(0) = 1, u'(0) = 0 → cos(t) + ode = Dt(Dt(u)) + u ~ 0 + @test ode isa GiacExpr + @test occursin("diff", string(ode)) + + u0 = u(0) ~ 1 + du0 = Dt(u)(0) ~ 0 + @test du0 isa DerivativeCondition + @test string(du0) == "u'(0)=0" + + result = desolve([ode, u0, du0], t, :u) + @test occursin("cos", string(result)) + end + + # ------------------------------------------------------------------ + # Symbolics parity (post-068 finishing touches) + # Convention: right-to-left Leibniz — in ∂ⁿf/∂v₁…∂vₙ the rightmost + # variable is applied first, the leftmost last. Consistent with + # operator product (∂/∂v₁)·…·(∂/∂vₙ) f. + # ------------------------------------------------------------------ + @testset "right-to-left convention in ∂-style show" begin + @giac_var x y z f(x, y) g(x, y, z) + + # Two-variable cross partial: Differential(y)(Differential(x)(f)) + # applies x first then y. Right-to-left convention places y + # (last-applied) leftmost in the denominator. + d_yx = Differential(y)(Differential(x)(f)) + io = IOBuffer(); show(io, d_yx) + @test String(take!(io)) == "D: ∂²f/∂y∂x" + + # Three-variable composition: applied x, y, z → denom is ∂z∂y∂x + d_zyx = Differential(z)(Differential(y)(Differential(x)(g))) + io = IOBuffer(); show(io, d_zyx) + @test String(take!(io)) == "D: ∂³g/∂z∂y∂x" + + # Higher order with two vars: steps = [("x", 2), ("y", 1)] + # applied x twice first, then y once → denom is ∂y∂x² + d_y_x2 = Differential(y)(Differential(x, 2)(f)) + @test d_y_x2.steps == [("x", 2), ("y", 1)] + io = IOBuffer(); show(io, d_y_x2) + @test String(take!(io)) == "D: ∂³f/∂y∂x²" + end + + @testset "Base.^(::Differential, ::Integer)" begin + @giac_var x t f(x) u(t) + + # Differential(t)^2 ≡ Differential(t, 2): same effect on a function-form + Dt2_pow = Differential(t)^2 + Dt2_ctor = Differential(t, 2) + @test Dt2_pow isa Differential + @test Dt2_pow.var === Dt2_ctor.var + @test Dt2_pow.order == Dt2_ctor.order == 2 + @test Dt2_pow(u).steps == Dt2_ctor(u).steps == [("t", 2)] + + # Higher exponent + @test (Differential(x)^3).order == 3 + @test (Differential(x)^3)(f).steps == [("x", 3)] + + # Composition: (Dx^2) * (Dx^1) chained as Dx^2 applied on a Dx step + # collapses thanks to _append_step, giving steps = [("x", 3)] + d = Differential(x)(f) + @test (Differential(x)^2)(d).steps == [("x", 3)] + + # n == 0 returns identity (Symbolics convention) + @test Differential(t)^0 === identity + @test (Differential(t)^0)(u) === u + + # n == 1 is a no-op on order + @test (Differential(t)^1).order == 1 + end + + @testset "Base.:*(::Differential, ::Differential) — composition" begin + @giac_var x y f(x, y) + + # (Dy * Dx)(f) ≡ Dy(Dx(f)) + composed = Differential(y) * Differential(x) + @test composed(f).steps == Differential(y)(Differential(x)(f)).steps == + [("x", 1), ("y", 1)] + + # Same variable: (Dx * Dx)(f) collapses to steps [("x", 2)] + @test (Differential(x) * Differential(x))(f).steps == [("x", 2)] + end + + @testset "Base.∘(::Differential, ::Differential) — composition" begin + @giac_var x y f(x, y) + + # (Dy ∘ Dx)(f) ≡ Dy(Dx(f)) + composed = Differential(y) ∘ Differential(x) + @test composed(f).steps == [("x", 1), ("y", 1)] + + # Mixing ∘ and * gives the same result + @test (Differential(y) ∘ Differential(x))(f).steps == + (Differential(y) * Differential(x))(f).steps + end + + @testset "expand_derivatives" begin + @giac_var x y t f(x, y) u(t) + + # On a plain GiacExpr → identity (already eager-evaluated in Giac.jl) + e = x^2 + y * x + @test expand_derivatives(e) === e + + # On the result of bare-expression Differential (already a GiacExpr) → identity + already = Differential(x)(x^2 + y * x) + @test already isa GiacExpr + @test expand_derivatives(already) === already + + # On a DerivativeExpr (function-form, lazy) → returns a plain GiacExpr. + # GIAC may reformulate `diff(f(x,y),x)` internally as `diff(f,0)(x,y)` + # (its operator-then-evaluate normal form) — semantically equivalent + # but syntactically different — so we only assert the type and that + # GIAC's diff was actually called (output mentions `diff` and `f`). + lazy = Differential(x)(f) + @test lazy isa DerivativeExpr + evaluated = expand_derivatives(lazy) + @test evaluated isa GiacExpr + @test !(evaluated isa DerivativeExpr) + s = string(evaluated) + @test occursin("diff", s) + @test occursin("f", s) + + # Differential(t, 2)(u) lazy → expand returns a GiacExpr referring to u + # (e.g. `diff(u(t),t,2)` or `diff(u,0,0)(t)` after GIAC normalization). + d2 = Differential(t, 2)(u) + @test d2 isa DerivativeExpr + ev2 = expand_derivatives(d2) + @test ev2 isa GiacExpr + @test occursin("diff", string(ev2)) + @test occursin("u", string(ev2)) + end + +end + diff --git a/test/test_types.jl b/test/test_types.jl index 060eba5c..269b0943 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -308,6 +308,26 @@ @test Giac._parse_function_expr("ψ(t)") == ("ψ", "t") @test Giac._parse_function_expr("α(β,γ)") == ("α", "β") end + + # spec 068: _parse_function_args returns full argument list + @testset "_parse_function_args returns full arg list (068-multivar-d-operator)" begin + @test Giac._parse_function_args("u(t)") == ("u", ["t"]) + @test Giac._parse_function_args("f(x,y)") == ("f", ["x", "y"]) + @test Giac._parse_function_args("g(a, b, c)") == ("g", ["a", "b", "c"]) + @test Giac._parse_function_args("h( x , y )") == ("h", ["x", "y"]) + + # Non-function expressions + @test Giac._parse_function_args("x") === nothing + @test Giac._parse_function_args("a+b") === nothing + + # GIAC operations are excluded + @test Giac._parse_function_args("diff(u,t)") === nothing + @test Giac._parse_function_args("sin(x)") === nothing + + # Unicode identifiers + @test Giac._parse_function_args("ϕ(𝑧)") == ("ϕ", ["𝑧"]) + @test Giac._parse_function_args("α(β,γ)") == ("α", ["β", "γ"]) + end end @testset "Derivative Operator D - Unicode identifiers" begin @@ -316,12 +336,11 @@ @giac_var 𝑧 ϕ(𝑧) dϕ = D(ϕ) @test dϕ isa DerivativeExpr - @test dϕ.order == 1 + @test dϕ.steps == [("𝑧", 1)] @test dϕ.funcname == "ϕ" - @test dϕ.varname == "𝑧" d2ϕ = D(ϕ, 2) - @test d2ϕ.order == 2 + @test d2ϕ.steps == [("𝑧", 2)] @test d2ϕ.funcname == "ϕ" end @@ -330,31 +349,30 @@ @giac_var u(t) du = D(u) @test du isa DerivativeExpr - @test du.order == 1 + @test du.steps == [("t", 1)] @test du.funcname == "u" - @test du.varname == "t" end @testset "D(u, 2) creates second derivative" begin @giac_var u(t) d2u = D(u, 2) @test d2u isa DerivativeExpr - @test d2u.order == 2 + @test d2u.steps == [("t", 2)] @test d2u.funcname == "u" end - @testset "D(D(u)) creates second derivative" begin + @testset "D(D(u)) creates second derivative (steps collapsed)" begin @giac_var u(t) d2u = D(D(u)) @test d2u isa DerivativeExpr - @test d2u.order == 2 + @test d2u.steps == [("t", 2)] end @testset "D(u, 3) creates third derivative" begin @giac_var y(t) d3y = D(y, 3) @test d3y isa DerivativeExpr - @test d3y.order == 3 + @test d3y.steps == [("t", 3)] end @testset "D on non-function raises error" begin @@ -475,5 +493,233 @@ show(io, D(u, 2)) @test String(take!(io)) == "D: u''(t)" end + + # spec 068 — n-th derivative notation: switch from primes to ⁽ⁿ⁾ for n ≥ 4 + @testset "show method uses ⁽ⁿ⁾ for n ≥ 4 (068-multivar-d-operator)" begin + @giac_var u(t) + + # Boundary: n = 3 still uses three primes + io = IOBuffer() + show(io, D(u, 3)) + @test String(take!(io)) == "D: u'''(t)" + + # Threshold: n = 4 switches to ⁽⁴⁾ notation + io = IOBuffer() + show(io, D(u, 4)) + @test String(take!(io)) == "D: u⁽⁴⁾(t)" + + # n = 5 + io = IOBuffer() + show(io, D(u, 5)) + @test String(take!(io)) == "D: u⁽⁵⁾(t)" + + # Two-digit superscript + io = IOBuffer() + show(io, D(u, 12)) + @test String(take!(io)) == "D: u⁽¹²⁾(t)" + end + + # spec 068 — DerivativePoint also uses the math-convention notation in show, + # but Base.string keeps prime notation (GIAC consumes prime form). + @testset "DerivativePoint show vs string (068-multivar-d-operator)" begin + @giac_var u(t) + + # n = 1 — show and string agree (single prime) + dp1 = D(u)(0) + @test string(dp1) == "u'(0)" + io = IOBuffer() + show(io, dp1) + @test String(take!(io)) == "u'(0)" + + # n = 4 — show uses ⁽⁴⁾, string keeps primes (for GIAC) + dp4 = D(u, 4)(0) + @test string(dp4) == "u''''(0)" # GIAC-consumable + io = IOBuffer() + show(io, dp4) + @test String(take!(io)) == "u⁽⁴⁾(0)" # human-readable + end + end + + # spec 068 US5 — transitional D(f, x) alias overloads + @testset "D transitional alias D(f, x) (068-multivar-d-operator US5)" begin + @giac_var x y f(x, y) + @giac_var t u(t) + + # D(f, x) ≡ Differential(x)(f), with depwarn + d1 = D(f, x) + d2 = Differential(x)(f) + @test d1.steps == d2.steps + @test d1.funcname == d2.funcname + + # D(f, x, 2) ≡ Differential(x)(Differential(x)(f)) + d3 = D(f, x, 2) + d4 = Differential(x)(Differential(x)(f)) + @test d3.steps == d4.steps + + # D(D(f, x), y) ≡ Differential(y)(Differential(x)(f)) + d5 = D(D(f, x), y) + d6 = Differential(y)(Differential(x)(f)) + @test d5.steps == d6.steps + + # D(u, t) ≡ Differential(t)(u) + d7 = D(u, t) + d8 = Differential(t)(u) + @test d7.steps == d8.steps + + # D(u, t, 3) ≡ third-order derivative + d9 = D(u, t, 3) + @test d9.steps == [("t", 3)] + + # Each transitional alias emits a depwarn (shape check on D(f, x)) + @test_logs (:warn, r"D\(expr, var\) is deprecated.*Differential") match_mode=:any D(f, x) + end + + # spec 068 US3 — D deprecation warnings + value-equivalence with Differential + @testset "D deprecation (068-multivar-d-operator US3)" begin + @giac_var t u(t) + + @testset "D(u) emits depwarn pointing to Differential(t)(u)" begin + # @test_logs captures logs (including Base.depwarn) regardless of --depwarn flag. + @test_logs (:warn, r"D\(u\) is deprecated.*Differential\(t\)\(u\)") match_mode=:any D(u) + end + + @testset "D(u, 2) emits depwarn referencing Differential" begin + @test_logs (:warn, r"D\(u, 2\) is deprecated.*Differential\(t\)") match_mode=:any D(u, 2) + end + + @testset "D(D(u)) emits depwarn on the chained call" begin + @test_logs (:warn, r"deprecated") match_mode=:any D(D(u)) + end + + @testset "deprecated D produces same value as canonical Differential" begin + # Result equality bypasses depwarn-shape checks: just compare structures. + d_via_D = D(u) + d_via_Differential = Differential(t)(u) + @test d_via_D.steps == d_via_Differential.steps + @test d_via_D.funcname == d_via_Differential.funcname + + d2_via_D = D(u, 2) + d2_via_Differential = Differential(t)(Differential(t)(u)) + @test d2_via_D.steps == d2_via_Differential.steps + + d2_chain_D = D(D(u)) + @test d2_chain_D.steps == d2_via_Differential.steps + end + + @testset "warn-once-per-callsite (Base.depwarn dedupe)" begin + # Base.depwarn deduplicates by (funcsym, file, line). Two calls from + # the same Julia line produce only one warning per session. We rely on + # Julia's documented Base.depwarn behavior; this test simply asserts + # that repeated D calls at the same site return correct values without + # error (warning emission shape covered by other tests above). + f() = (D(u); D(u)) + d1, d2 = f(), f() + @test d1 isa DerivativeExpr + @test d2 isa DerivativeExpr + @test d1.steps == d2.steps == [("t", 1)] + end + end + + # spec 068 US2 — multi-arg D(f) ambiguity guard + @testset "D(f) ambiguity guard on multi-arg function (068-multivar-d-operator US2)" begin + @giac_var x y f(x, y) + # D(f) on multi-arg function must raise ArgumentError, not silently differentiate. + err = nothing + try + D(f) + catch e + err = e + end + @test err isa ArgumentError + msg = err.msg + @test occursin("f", msg) + @test occursin("x", msg) + @test occursin("y", msg) + @test occursin("Differential(x)(f)", msg) + @test occursin("deprecated", msg) + + # D(f, n) on multi-arg also raises with the same shape + err2 = nothing + try + D(f, 2) + catch e + err2 = e + end + @test err2 isa ArgumentError + @test occursin("multiple arguments", err2.msg) + @test occursin("Differential", err2.msg) + + # Three-variable case + @giac_var a b c g(a, b, c) + err3 = nothing + try + D(g) + catch e + err3 = e + end + @test err3 isa ArgumentError + msg3 = err3.msg + @test occursin("g", msg3) + @test occursin("a", msg3) + @test occursin("b", msg3) + @test occursin("c", msg3) + + # Mono-var D(u) still succeeds (regression guard) + @giac_var u(t) + @test D(u) isa DerivativeExpr + end + + # spec 068 — foundational refactor of DerivativeExpr (multi-step) + @testset "DerivativeExpr multi-step (068-multivar-d-operator)" begin + @giac_var x y f(x, y) + + @testset "manual multi-step construction → string is nested diff" begin + d_xy = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 1), ("y", 1)]) + @test string(d_xy) == "diff(diff(f(x,y),x),y)" + + d_x2y = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 2), ("y", 1)]) + @test string(d_x2y) == "diff(diff(f(x,y),x,2),y)" + end + + @testset "mono-step DerivativeExpr keeps single diff form" begin + @giac_var u(t) + @test string(D(u)) == "diff(u(t),t)" + @test string(D(u, 2)) == "diff(u(t),t,2)" + end + + @testset "multi-step show uses ∂-style notation (right-to-left)" begin + # Convention: in ∂ⁿf/∂v₁…∂vₙ the rightmost variable is applied first, + # the leftmost last — consistent with the operator product + # (∂/∂v₁)·…·(∂/∂vₙ) f. Therefore `steps` (innermost-first) is + # reversed in the denominator. + d_xy = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 1), ("y", 1)]) + io = IOBuffer() + show(io, d_xy) + @test String(take!(io)) == "D: ∂²f/∂y∂x" + + d_x2y = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 2), ("y", 1)]) + io = IOBuffer() + show(io, d_x2y) + @test String(take!(io)) == "D: ∂³f/∂y∂x²" + end + + @testset "(d::DerivativeExpr)(args...) raises on multi-step" begin + d_xy = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 1), ("y", 1)]) + @test_throws ArgumentError d_xy(0, 0) + end + + @testset "D(::DerivativeExpr) raises on multi-step" begin + d_xy = DerivativeExpr(f, "f", Tuple{String, Int}[("x", 1), ("y", 1)]) + @test_throws ArgumentError D(d_xy) + @test_throws ArgumentError D(d_xy, 2) + end + + @testset "_to_giac_expr on mono-step still works (arithmetic surface)" begin + @giac_var u(t) + # If this round-trips through giac_eval correctly, arithmetic ops keep working. + result = D(u) + u + @test result isa GiacExpr + @test occursin("diff", string(result)) + end end end