56 Symbolics.jl

There are a few options in Julia for symbolic math, for example, the SymPy package which wraps a Python library. This section describes a collection of native Julia packages providing many features of symbolic math.

56.1 About

The Symbolics package bills itself as a “fast and modern Computer Algebra System (CAS) for a fast and modern programming language.” This package relies on the SymbolicUtils package and is built upon by the ModelingToolkit package, which we don’t describe here.

We begin by loading the Symbolics package which when loaded re-exports the SymbolicUtils package.

using Symbolics

56.2 Symbolic variables

Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two packages discussed here.

SymbolicUtils creates variables which carry Julia type information (e.g. Int, Float64, …). This type information carries through operations involving these variables. Symbolic variables can be created with the @syms macro. For example:

@syms x y::Int f(x::Real)::Real

(x, y, f)

This creates x a symbolic value with symbolic type Number, y a symbolic variable holding integer values, and f a symbolic function of a single real variable outputting a real variable.

The non-exported symtype function reveals the underlying type:

import Symbolics.SymbolicUtils: symtype

symtype(x), symtype(y)

(Number, Int64)

For y, the symbolic type being real does not imply the type of y is a subtype of Real:

isa(y, Real)

false

We see that the function f when called with y would return a value of (symbolic) type Real:

f(y) |> symtype

Real

As the symbolic type of x is Number – which is not a subtype of Real – the following will error:

f(x)

Tuple{Number} is not a subtype of Tuple{Real}.
Stacktrace:
 [1] error(s::LazyString)
   @ Base ./error.jl:44
 [2] promote_symtype(f::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, args::Type)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/types.jl:1080
 [3] (::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})(args::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/types.jl:1006
 [4] top-level scope
   @ ~/julia/CalculusWithJuliaNotes/quarto/alternatives/symbolics.qmd:64

Note

The SymPy package also has an @syms macro to create variables. Though their names agree, they do different things. Using both packages together would require qualifying many shared method names. For SymbolicUtils, the @syms macro uses Julia types to parameterize the variables. In SymPy it is possible to specify assumptions on the variables, but that is different and not useful for dispatch without some extra effort.

For Symbolics, symbolic variables are created using a wrapper around an underlying SymbolicUtils object. This wrapper, Num, is a subtype of Real (the underlying SymbolicUtils object may have symbolic type Real, but it won’t be a subtype of Real.)

Symbolic values are created with the @variables macro. For example:

@variables x y::Int z[1:3]::Int f(..)::Int

4-element Vector{Any}:
 x
 y
  z[1:3]
  f⋆

This creates

a symbolic value x of symtype Real
a symbolic value y of symtype Int
a vector of symbolic values each of symtype Int
a symbolic function f returning an object of symtype Int

The symbolic type reflects that of the underlying object behind the Num wrapper:

typeof(x), symtype(x), typeof(Symbolics.value(x))

(Num, Real, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})

(The value method unwraps the Num wrapper.)

56.3 Symbolic expressions

Symbolic expressions are built up from symbolic variables through natural Julia idioms. SymbolicUtils privileges a few key operations: Add, Mul, Pow, and Div. For example:

@syms x y
typeof(x + y) # `Add`

SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}

typeof(x * y) # `Mul`

SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}

Whereas, applying a function leaves a different type:

typeof(sin(x))

SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}

The Term wrapper just represents the effect of calling a function (in this case sin) on its arguments (in this case x).

This happens in the background with symbolic variables in Symbolics:

@variables x
typeof(sin(x)), typeof(Symbolics.value(sin(x)))

(Num, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})

56.3.1 Tree structure to expressions

The TermInterface package is used by SymbolicUtils to explore the tree structure of an expression. The main methods are (cf. SymbolicUtils.jl):

iscall(ex): true if ex is not a leaf node (like a symbol or numeric literal). The old name was istree.
operation(ex): the function being called (if iscall returns true)
arguments(ex): the arguments to the function being called
symtype(ex): the inferred type of the expression

In addition, the issym function, to determine if x is of type Sym, is useful to distinguish leaf nodes, as will be illustrated below.

These methods can be used to “walk” the tree:

@syms x y
ex = 1 + x^2 + y
operation(ex) # the outer function is `+`

+ (generic function with 288 methods)

arguments(ex) # `+` is n-ary, in this case with 3 arguments

3-element ReadOnlyArrays.ReadOnlyVector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.SmallVec{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}}:
 1
 x^2
 y

ex1 = arguments(ex)[2] # terms have been reordered
operation(ex1)  # operation for `x^2` is `^`

^ (generic function with 100 methods)

a, b = arguments(ex1)

2-element ReadOnlyArrays.ReadOnlyVector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.SmallVec{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}}:
 x
 2

iscall(ex1), iscall(a)

(true, false)

Here a is not a call, as it has no operation or arguments, it is just a variable (the x variable).

The value of symtype is the inferred type of an expression, which may not match the actual type. For example,

@variables x::Int
symtype(x), symtype(sin(x)), symtype(x/x), symtype(x / x^2)

(Int64, Real, Int64, Real)

The last one, is not likely to be an integer, but that is the inferred type in this case.

Example

As an example, we write a function to find the free symbols in a symbolic expression comprised of SymbolicUtils variables. (The Symbolics.get_variables also does this task.) To find the symbols involves walking the expression tree until a leaf node is found and then adding that to our collection if it matches issym.

import Symbolics.SymbolicUtils: issym
free_symbols(ex) = (s=Set(); free_symbols!(s, ex); s)
function free_symbols!(s, ex)
    if iscall(ex)
        for a ∈ arguments(ex)
            free_symbols!(s, a)
        end
    else
        issym(ex) && push!(s, ex) # push new symbol onto set
    end
end

free_symbols! (generic function with 1 method)

@syms x y z
ex = sin(x + 1)*cos(z)
free_symbols(ex)

Set{Any} with 2 elements:
  z
  x

56.4 Expression manipulation

56.4.1 Substitute

The substitute command is used to replace values with other values. For example:

@variables x y z
ex = 1 + x + x^2/2 + x^3/6
substitute(ex, x=>1)

8//3

This defines a symbolic expression, then substitutes the value 1 in for x. The Pair notation is useful for a single substitution. When there is more than one substitution, a dictionary is used:

w = x^3 + y^3 - 2z^3
substitute(w, Dict(x=>2, y=>3))

35 - 2(z^3)

56.4.2 Simplify

Algebraic operations with symbolic values can involve an exponentially increasing number of terms. As such, some simplification rules are applied after an operation to reduce the complexity of the computed value.

For example, 0+x should simplify to x, as well 1*x, x^0, or x^1 should each simplify, to some natural answer.

SymbolicUtils also simplifies several other expressions, including:

-x becomes (-1)*x
x * x becomes x^2 (and x^n if more terms). Meaning this expression is represented as a power, not a product
x + x becomes 2*x (and n*x if more terms). Similarly, this is represented as a product, not a sum.
p/q * x becomes (p*x)/q), similarly p/q * x/y becomes (p*x)/(q*y). (Division wraps multiplication.)

In SymbolicUtils, this rewriting is accomplished by means of rewrite rules. The package makes it easy to apply user-written rewrite rules.

56.4.3 Rewriting

Many algebraic simplifications are done by the simplify command. For example, the basic trigonometric identities are applied:

@variables x
ex = sin(x)^2 + cos(x)^2
ex, simplify(ex)

(sin(x)^2 + cos(x)^2, 1)

The simplify function applies a series of rewriting rule until the expression stabilizes. The rewrite rules can be user generated, if desired. For example, the Pythagorean identity of trigonometry, just used, can be implemented with this rule:

r = @acrule(sin(~x)^2 + cos(~x)^2 => one(~x))
ex |> Symbolics.value |> r |> Num

The rewrite rule, r, is defined by the @acrule macro. The a is for associative, the c for commutative, assumptions made by the macro. (The c means cos(x)^2 + sin(x)^2 will also simplify.) Rewrite rules are called on the underlying SymbolicUtils expression, so we first unwrap, then after re-wrap.

The above expression for r is fairly easy to appreciate. The value ~x matches the same variable or expression. So the above rule will also simplify more complicated expressions:

@variables y z
ex1 = substitute(ex, x => sin(x + y + z))
ex1 |> Symbolics.value |> r |> Num

Rewrite rules when applied return the rewritten expression, if there is a match, or nothing.

Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trigonometry and exponential simplifications:

r = @rule(exp(~x)^(~y) => exp(~x * ~y))  # (e^x)^y -> e^(x*y)
ex = exp(-x+z)^y
ex, ex |> Symbolics.value |> r |> Num

(exp(-x + z)^y, exp((-x + z)*y))

This rule is not commutative or associative, as x^y is not the same as y^x and (x^y)^z is not x^(y^z) in general.

The application of rules can be filtered through qualifying predicates. This artificial example uses iseven which returns true for even numbers. Here we subtract 1 when a number is not even, and otherwise leave the number alone. We do this with two rules:

reven = @rule ~x::iseven => ~x
rodd = @rule ~x::(!iseven) => ~x - 1
r = SymbolicUtils.Chain([rodd, reven])
r(2), r(3)

(2, 2)

The Chain function conveniently allows the sequential application of rewrite rules.

The notation ~x is called a “slot variable” in the documentation for SymbolicUtils. It matches a single expression. To match more than one expression, a “segment variable”, denoted with two ~s is used.

56.4.4 Creating functions

By utilizing the tree-like nature of a symbolic expression, a Julia expression can be built from a symbolic expression easily enough. The Symbolics.toexpr function does this:

ex = exp(-x + z)^y
Symbolics.toexpr(ex)

:((^)((exp)((+)((*)(-1, x), z)), y))

This output shows an internal representation of the steps for computing the value ex given different inputs. (The number (-1) multiplies x, this is added to z and the result passed to exp. That values is then used as the base for ^ with exponent y.)

Such Julia expressions are one step away from building Julia functions for evaluating symbolic expressions fast (though with some technical details about “world age” to be reckoned with). The build_function function with the argument expression=Val(false) will compile a Julia function:

h = build_function(ex, x, y, z; expression=Val(false))
h(1, 2, 3)

54.59815003314424

The above is similar to substitution:

substitute(ex, Dict(x=>1, y=>2, z=>3))

exp(2)^2

However, build_function will be significantly more performant, which when many function calls are used – such as with plotting – is a big advantage.

Note

The documentation colorfully says “build_function is kind of like if lambdify (from SymPy) ate its spinach.”

The above, through passing \(3\) variables after the expression, creates a function of \(3\) variables. Functions of a vector of inputs can also be created, just by expressing the variables in that manner:

h1 = build_function(ex, [x, y, z]; expression=Val(false))
h1([1, 2, 3])  # not h1(1,2,3)

54.59815003314424

Example

As an example, here we use the Roots package to find a zero of a function defined symbolically:

import Roots
@variables x
ex = x^5 - x - 1
λ = build_function(ex, x; expression=Val(false))
Roots.find_zero(λ, (1, 2))

1.1673039782614187

56.5 Plotting

Using Plots, the plotting of symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via build_function.

For example,

using Plots
@variables x
plot(x^x^x, 0, 2)

A parametric plot is easily defined:

plot(sin(x), cos(x), 0, pi/4)

Expressions to be plotted can represent multivariate functions.

@variables x y
ex = 3*(1-x)^2*exp(-x^2 - (y+1)^2) - 10(x/5-x^3-y^5)*exp(-x^2-y^2) - 1/3*exp(-(x+1)^2-y^2)
xs = ys = range(-5, 5, length=100)
surface(xs, ys, ex)

The ordering of the variables is determined by Symbolics.get_variables:

Symbolics.get_variables(ex)

OrderedCollections.OrderedSet{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}} with 2 elements:
  x
  y

56.6 Polynomial manipulations

There are some facilities for manipulating polynomial expressions in Symbolics. A polynomial, mathematically, is an expression involving one or more symbols with coefficients from a collection that has, at a minimum, addition and multiplication defined. The basic building blocks of polynomials are monomials, which are comprised of products of powers of the symbols. Mathematically, monomials are often allowed to have a multiplying coefficient and may be just a coefficient (if each symbol is taken to the power \(0\)), but here we consider just expressions of the type \(x_1^{a_1} \cdot x_2^{a_2} \cdots \cdot x_k^{a_k}\) with the \(a_i > 0\) as monomials.

With this understanding, then an expression can be broken up into monomials with a possible coefficient (possibly just \(1\)) and terms which are not monomials (such as a constant or a more complicated function of the symbols). This is what is returned by the polynomial_coeffs function.

For example

@variables a b c x
d, r = polynomial_coeffs(a*x^2 + b*x + c, (x,))

(Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}(1 => c, x^2 => a, x => b), 0)

The first term output is a dictionary with keys which are the monomials and with values which are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just \(0\).

The expression can then be reconstructed through

r + sum(v*k for (k,v) ∈ d)

c + b*x + a*(x^2)

The above has a,b,c as parameters and x as the symbol. This separation is designated by passing the desired polynomial symbols to polynomial_coeff as an iterable. (Above as a \(1\)-element tuple.)

More complicated polynomials can be similarly decomposed:

@variables a b c x y z
ex = a*x^2*y*z + b*x*y^2*z + c*x*y*z^2
d, r = polynomial_coeffs(ex, (x, y, z))

(Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}(x*(y^2)*z => b, (x^2)*y*z => a, x*y*(z^2) => c), 0)

The (sparse) decomposition of the polynomial is returned through d. The same pattern as above can be used to reconstruct the expression. To extract the coefficient for a monomial term, indexing can be used. Of note, is an expression like x^2*y*z could possibly not equal the algebraically equal x*y*z*x, as they are only equal after some simplification, but the keys are in a canonical form, so this is not a concern:

d[x*y*z*x], d[z*y*x^2]

(a, a)

The residual term will capture any non-polynomial terms:

ex = sin(x) - x + x^3/6
d, r = polynomial_coeffs(ex, (x,))
r

sin(x)

To find the degree of a monomial expression, the degree function is available, though not exported. Here it is applied to each monomial in d:

import Symbolics: degree
[degree(k) for (k,v) ∈ d]

2-element Vector{Int64}:
 1
 3

The degree function will also identify the degree of more complicated terms:

degree(1 + x + x^2)

Mathematically the degree of the \(0\) polynomial may be \(-1\) or undefined, but here it is \(0\):

degree(0), degree(1), degree(x), degree(x^a)

(0, 0, 1, a)

The coefficients are returned as values of a dictionary, and dictionaries are unsorted.

@variables x a0 as[1:10]
p = a0 + sum(as[i]*x^i for i ∈ eachindex(collect(as)))
d, r = polynomial_coeffs(p, (x,))
d

Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}} with 11 entries:
  x^10 => as[10]
  x    => as[1]
  1    => a0
  x^7  => as[7]
  x^5  => as[5]
  x^3  => as[3]
  x^6  => as[6]
  x^4  => as[4]
  x^8  => as[8]
  x^2  => as[2]
  x^9  => as[9]

To have a natural map between polynomials of a single symbol in the standard basis and a vector, we could use a pattern like this to sort the values:

vcat(r, [d[k] for k ∈ sort(collect(keys(d)), by=degree)])

12-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 0
 a0
 as[1]
 as[2]
 as[3]
 as[4]
 as[5]
 as[6]
 as[7]
 as[8]
 as[9]
 as[10]

As an example usage, we write a function that can determine if an expression is a polynomial expression over some specified variables.

function is_poly(expr, vars)
    all(Symbolics.SymbolicUtils.issym ∘ Symbolics.value, vars) || error("vars must be an iterable of symbols")
    p,r = polynomial_coeffs(expr, vars)
    length(intersect(Symbolics.get_variables(r), vars)) == 0
end

is_poly(p, (x,))

true

Rational expressions can be decomposed into a numerator and denominator using the following idiom, which assumes the outer operation is division (a binary operation):

@variables x
ex = (1 + x + x^2) / (1 + x + x^2 + x^3)
function nd(ex)
    ex1 = Symbolics.value(ex)
    (operation(ex1) == /) || return (ex, one(ex))
    Num.(arguments(ex1))
end
nd(ex)

2-element Vector{Num}:
       1 + x + x^2
 1 + x + x^2 + x^3

With this, the study of asymptotic behaviour of a univariate rational expression would involve an investigation like the following:

m,n = degree.(nd(ex))
m > n ? "limit is infinite" : m < n ? "limit is 0" : "limit is a constant"

"limit is 0"

56.7 Vectors and matrices

Symbolic vectors and matrices can be created with a specified size:

@variables v[1:3] M[1:2, 1:3] N[1:3, 1:3]

3-element Vector{Symbolics.Arr{Num}}:
 v[1:3]
 M[1:2,1:3]
 N[1:3,1:3]

Computations, like finding the determinant below, are lazy unless the values are collected:

using LinearAlgebra
det(N)

det(N)

det(collect(N))

N[1, 1]*(N[2, 2]*N[3, 3] - N[2, 3]*N[3, 2]) - N[1, 2]*(N[2, 1]*N[3, 3] - N[2, 3]*N[3, 1]) + N[1, 3]*(N[2, 1]*N[3, 2] - N[2, 2]*N[3, 1])

Similarly, with norm, which returns a generator unless collected:

norm(collect(v))

sqrt(abs2(v[1]) + abs2(v[2]) + abs2(v[3]))

Matrix multiplication is also deferred, but the size compatibility of the matrices and vectors is considered immediately:

M*N, N*N, M*v

(M*N, N^2, M*v)

This errors, as the matrix dimensions are not compatible for multiplication:

N*M

ArgumentError: Encountered incompatible shapes UnitRange{Int64}[1:3, 1:3] and UnitRange{Int64}[1:2, 1:3] when multiplying.
Stacktrace:
 [1] throw_incompatible_shapes(x::SymbolicUtils.SmallVec{UnitRange{Int64}, Vector{UnitRange{Int64}}}, y::SymbolicUtils.SmallVec{UnitRange{Int64}, Vector{UnitRange{Int64}}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:10
 [2] _multiplied_shape(shapes::Tuple{SymbolicUtils.SmallVec{UnitRange{Int64}, Vector{UnitRange{Int64}}}, SymbolicUtils.SmallVec{UnitRange{Int64}, Vector{UnitRange{Int64}}}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:47
 [3] _multiplied_terms_shape(terms::Tuple{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:83
 [4] (::SymbolicUtils.MulWorkerBuffer{SymReal})(terms::Tuple{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:217
 [5] mul_worker(::Type{SymReal}, terms::Tuple{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:357
 [6] *
   @ ~/.julia/packages/SymbolicUtils/XxhMQ/src/symbolic_ops/mul.jl:361 [inlined]
 [7] *(::Symbolics.Arr{Num, 2}, ::Symbolics.Arr{Num, 2})
   @ Symbolics ~/.julia/packages/Symbolics/7hRc0/src/array-lib.jl:103
 [8] top-level scope
   @ ~/julia/CalculusWithJuliaNotes/quarto/alternatives/symbolics.qmd:600

Similarly, linear solutions can be symbolically specified:

@variables R[1:2, 1:2] b[1:2]
R \ b

(R \ b)[1:2]

collect(R \ b)

2-element Vector{Num}:
 (R \ b)[1]
 (R \ b)[2]

56.8 Algebraically solving equations

The ~ operator creates a symbolic equation. For example

@variables x y
x^5 - x ~ 1

-x + x^5 ~ 1

eqs = [5x + 2y, 6x + 3y] .~ [1, 2]

2-element Vector{Equation}:
 5x + 2y ~ 1
 6x + 3y ~ 2

The Symbolics.symbolic_linear_solve function can solve linear equations. For example,

Symbolics.symbolic_linear_solve(eqs, [x, y])

2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 (-1//3)
 (4//3)

The coefficients can be symbolic. Two examples could be:

@variables m b x y
eq = y ~ m*x + b
Symbolics.symbolic_linear_solve(eq, x)

(-b + y) / m

@variables a11 a12 a22 x y b1 b2
R,X,b = [a11 a12; 0 a22], [x; y], [b1, b2]
eqs = R*X .~ b

2-element Vector{Equation}:
 a11*x + a12*y ~ b1
 a22*y ~ b2

Symbolics.symbolic_linear_solve(eqs, [x,y])

2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}:
 (-b1 + (a12*b2) / a22) / (-a11)
 b2 / a22

56.9 Limits

Many symbolic limits involving exponentials and logarithms can be computed in Symbolics, as of recent versions. The underlying package is SymbolicLimits. This package is in development. Below we use the unwrapped version of the variable. We express limits as \(x\) goes to infinity, which can be achieved by rewriting:

@variables x
𝑥 = x.val  # unwrapped
F(x) = exp(x+exp(-x))-exp(x)
limit(F(𝑥), 𝑥, Inf)

1.0

F(x) = log(log(x*exp(x*exp(x))+1))-exp(exp(log(log(x))+1/x))
limit(F(𝑥), 𝑥, Inf)

56.10 Derivatives

Symbolics provides the derivative function to compute the derivative of a function with respect to a variable:

@variables a b c x
y = a*x^2 + b*x + c
yp = Symbolics.derivative(y, x)

b + 2a*x

Or to find a critical point:

Symbolics.symbolic_linear_solve(yp ~ 0, x) # linear equation to solve

b / (-2a)

The derivative computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:

D = Differential(x)
D(y)

Differential(x, 1)(c + b*x + a*(x^2))

and then

expand_derivatives(D(y))

b + 2a*x

Using Differential, differential equations can be specified. An example was given in ODEs, using ModelingToolkit.

Higher order derivatives can be done through composition:

D(D(y)) |> expand_derivatives

2a

Differentials can also be multiplied to create operators for taking higher-order derivatives:

@variables x y
ex = (x - y^2)/(x^2 + y^2)
Dx, Dy = Differential(x), Differential(y)
Dxx, Dxy, Dyy = Dx*Dx, Dx*Dy, Dy*Dy
[Dxx(ex) Dxy(ex); Dxy(ex) Dyy(ex)] .|> expand_derivatives

2×2 Matrix{Num}:
 (8(x^2)*(x - (y^2))) / ((x^2 + y^2)^3) + (-2x) / ((x^2 + y^2)^2) + (-2x - 2(x - (y^2))) / ((x^2 + y^2)^2)  …                                           (4x*y) / ((x^2 + y^2)^2) + (-2y) / ((x^2 + y^2)^2) + (8x*y*(x - (y^2))) / ((x^2 + y^2)^3)
                 (4x*y) / ((x^2 + y^2)^2) + (-2y) / ((x^2 + y^2)^2) + (8x*y*(x - (y^2))) / ((x^2 + y^2)^3)     (8(x - (y^2))*(y^2)) / ((x^2 + y^2)^3) + (4(y^2)) / ((x^2 + y^2)^2) + (4(y^2) - 2(x - (y^2))) / ((x^2 + y^2)^2) + -2 / (x^2 + y^2)

In addition to Symbolics.derivative there are also the helper functions, such as hessian which performs the above

Symbolics.hessian(ex, [x,y])

2×2 Matrix{Num}:
 (8(x^2)*(x - (y^2))) / ((x^2 + y^2)^3) + (-2x) / ((x^2 + y^2)^2) + (-2x - 2(x - (y^2))) / ((x^2 + y^2)^2)  …                                           (4x*y) / ((x^2 + y^2)^2) + (-2y) / ((x^2 + y^2)^2) + (8x*y*(x - (y^2))) / ((x^2 + y^2)^3)
                 (4x*y) / ((x^2 + y^2)^2) + (-2y) / ((x^2 + y^2)^2) + (8x*y*(x - (y^2))) / ((x^2 + y^2)^3)     (8(x - (y^2))*(y^2)) / ((x^2 + y^2)^3) + (4(y^2)) / ((x^2 + y^2)^2) + (4(y^2) - 2(x - (y^2))) / ((x^2 + y^2)^2) + -2 / (x^2 + y^2)

The gradient function is also defined

@variables x y z
ex = x^2 - 2x*y + z*y
Symbolics.gradient(ex, [x, y, z])

3-element Vector{Num}:
 2x - 2y
 -2x + z
       y

The jacobian function takes an array of expressions:

@variables x y
eqs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(eqs, [x,y])

2×2 Matrix{Num}:
 2x  -2y
 2y   2x

56.11 Integration

The SymbolicIntegration package provides two means to integration univariate functions using either the Risch algorithm or a rules-based approach.

using SymbolicIntegration, Symbolics

The main entry point is the function integrate.

This is a substitution test:

@variables x a b
integrate(x*exp(-x^2), x)

(ℯ^(-(x^2))) / -2

This is an integration by parts example:

integrate(x * sin(x), x)

sin(x) - x*cos(x)

As is this. We

integrate(a * log(b*x), x)

a*(-x + x*log(b*x))

The integration of rational functions (ratios of polynomials) can be done algorithmically, provided the underlying factorizations can be identified.

integrate((x-1)/((x-2)^3*(x-4)), x)

(3//8)*log(4 - x) + (-(5//4) + (3//4)*x) / (4 - 4x + x^2) - (3//8)*log(2 - x)

Method selection

A third argument to integrate can specify the method, as in integrate(expr, var, method).

This function does not get solved by the Risch method:

ex = cos(5x)*sin(x)
integrate(ex, x, RischMethod())

∫(cos(5x)*sin(x), x)

However, it can be done with the rules-based method:

integrate(ex, x, RuleBasedMethod())

-(1//12)*cos(6x) + (1//8)*cos(-4x)

The example of integrating (x-1)/((x-2)^3*(x-4)), done above, is one where the RischMethod works, but not the RuleBasedMethod.

Each method has different keyword arguments. For RuleBasedMethod the verbose=true argument will show which rules were applied. In this case, there is a single one:

integrate(ex, x, RuleBasedMethod(verbose=true))

┌-------Applied rule 4_7_4_3 on ∫(cos(5x)*sin(x), x)
| ∫ sin(a + b * x) * cos(c + d * x) dx => if 
|       !(contains_var(a, b, c, d, x)) &&
|       b ^ 2 - d ^ 2 !== 0
| -(cos((a - c) + (b - d) * x)) / (2 * (b - d)) - cos(a + c + (b + d) * x) / (2 * (b + d))
└-------with result: -(1//12)*cos(6x) + (1//8)*cos(-4x)

-(1//12)*cos(6x) + (1//8)*cos(-4x)

This example, shows two rules are applied:

integrate(a * log(b*x), x, RuleBasedMethod(verbose=true))

┌-------Applied rule 0_1_12 on ∫(a*log(b*x), x)
| ∫ a*f(x) dx => a*∫ f(x) dx
└-------with result: a*∫(log(b*x), x)
┌-------Applied rule 3_1_1_1 on ∫(log(b*x), x)
| ∫ log(c * x ^ n) dx => if 
|       !(contains_var(c, n, x))
| x * log(c * x ^ n) - n * x
└-------with result: -x + x*log(b*x)

a*(-x + x*log(b*x))

To read the rules, there are different predicate functions involved. The commonly used contains_var predicate checks through pattern matching if the last variable is contained in any of the rest of the specified variables. In the pattern c * x ^ n the match is c is b and n is 1. Neither depends on x, so this rule is applied to integrate log(b*x).

The contains_var(c, n, x) call above checks if the constant b (matching c in the pattern) depends on x and if 1 depends on x (matching the power n, using a default for the variable, in the pattern). The rules comes from Rubi which has some 7000 rules, many of which are implemented in SymbolicIntegration.