Equation solving

[ChrisRackauckas]

This looks like a good spot to implement equation solving. Similar to https://docs.sympy.org/latest/modules/solvers/solvers.html .

• [ ] Polynomials
• [ ] Nonlinear 1D equations
• [ ] Systems of equations

[dlfivefifty]

• [ ] Eigenvalues

julia> @syms x::Real
(x,)

julia> [0 x; x 0]
2×2 Array{Any,2}:
 0    x
  x  0

julia> [0 x; x 0] |> eigvals
ERROR: MethodError: no method matching zero(::Type{Any})
Closest candidates are:
  zero(::Type{Union{Missing, T}}) where T at missing.jl:105
  zero(::Type{Missing}) at missing.jl:103
  zero(::Type{LibGit2.GitHash}) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LibGit2/src/oid.jl:220
  ...
Stacktrace:
 [1] zero(::Type{Any}) at ./missing.jl:105
 [2] eigtype(::Type{T} where T) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:302
 [3] eigvals(::Array{Any,2}; kws::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326
 [4] eigvals(::Array{Any,2}) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326
 [5] |>(::Array{Any,2}, ::typeof(eigvals)) at ./operators.jl:823
 [6] top-level scope at REPL[8]:1
 [7] eval(::Module, ::Any) at ./boot.jl:331
 [8] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [9] run_backend(::REPL.REPLBackend) at /Users/solver/.julia/packages/Revise/MgvIv/src/Revise.jl:1023
 [10] top-level scope at none:0

[ChrisRackauckas]

Eigenvalues would need a polynomial equation solver. The generic fallbacks likely won't work

[dlfivefifty]

Mathematica's solution is a Root type to represent polynomial roots:

In[20]:= ( {
   {0, x, y, x},
   {x, 0, x, x},
   {0, x, x, x},
   {x, x, x, x}
  } ) // Eigenvalues

Out[20]= {Root[
  x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 2 x #1^3 + #1^4 &, 
  1], Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
    2 x #1^3 + #1^4 &, 2], 
 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
    2 x #1^3 + #1^4 &, 3], 
 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
    2 x #1^3 + #1^4 &, 4]}

[owiecc]

Is there a workaround to solve a system of nonlinear equations in a different engine and get the solution back to Julia? I already have a big set of equations modelled with Symbolics.jl and I just saw systems of nonlinear equations are not implemented yet. I'd prefer not to reimplement everything from scratch.

[anandijain]

Is there a workaround to solve a system of nonlinear equations in a different engine and get the solution back to Julia? I already have a big set of equations modelled with Symbolics.jl and I just saw systems of nonlinear equations are not implemented yet. I'd prefer not to reimplement everything from scratch.

@owiecc Symbolics.symbolics_to_sympy might be of help, not sure though

[YingboMa]

We could construct the characteristic polynomial and introduce the Root object

julia> using LinearAlgebra, Symbolics

julia> @variables x y λ
3-element Vector{Num}:
 x
 y
 λ

julia> A = [0 x y x; x 0 x x; 0 x x x; x x x x]
4×4 Matrix{Num}:
 0  x  y  x
 x  0  x  x
 0  x  x  x
 x  x  x  x

julia> simplify(det(A - I*λ), expand=true)
x^4 + λ^4 + λ*(x^3) - (2x*(λ^3)) - (y*(x^3)) - (4(x^2)*(λ^2)) - (2y*λ*(x^2))

Mathematica's root object is pretty smart. D is also defined on them.

In[10]:= ({{0, x, y, x}, {x, 0, x, x}, {0, x, x, x}, {x, x, x, x}}) //
   Eigenvalues;
D[%, x]

Out[11]= {(-4 x^3 + 
    3 x^2 y - (3 x^2 - 4 x y) Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1] + 
    8 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1]^2 + 
    2 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1]^3)/(x^3 - 2 x^2 y - 
    8 x^2 Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1] - 
    6 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1]^2 + 
    4 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 1]^3), (-4 x^3 + 
    3 x^2 y - (3 x^2 - 4 x y) Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2] + 
    8 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2]^2 + 
    2 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2]^3)/(x^3 - 2 x^2 y - 
    8 x^2 Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2] - 
    6 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2]^2 + 
    4 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 2]^3), (-4 x^3 + 
    3 x^2 y - (3 x^2 - 4 x y) Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3] + 
    8 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3]^2 + 
    2 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3]^3)/(x^3 - 2 x^2 y - 
    8 x^2 Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3] - 
    6 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3]^2 + 
    4 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 3]^3), (-4 x^3 + 
    3 x^2 y - (3 x^2 - 4 x y) Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4] + 
    8 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4]^2 + 
    2 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4]^3)/(x^3 - 2 x^2 y - 
    8 x^2 Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4] - 
    6 x Root[
      x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4]^2 + 
    4 Root[x^4 - x^3 y + (x^3 - 2 x^2 y) #1 - 4 x^2 #1^2 - 
        2 x #1^3 + #1^4 &, 4]^3)}

[blegat]

For systems of polynomials equations, there is an interface in SemialgebraicSets.jl. Two algorithms currently implement this interface:

• Compute Groebner basis, get the multiplication matrices and simultaneously diagonalize them with Schur factorization of a random linear combination of them. This is implemented in SemialgebraicSets.jl as well as the Groebner basis computation using Buchberger's algorithm. In the future, I would like to have other Groebner basis computation algorithm, e.g., by having https://github.com/ederc/GroebnerBasis.jl implement the interface.
• Homotopy continuation: https://www.juliahomotopycontinuation.org/

SymbolicUtils already has some code to transform symbolic expressions into multivariate polynomials using MultivariatePolynomials.jl so it should be too hard to link to SemialgebraicSets.jl.

[ChrisRackauckas]

GroebnerBasis.jl is a problematic dependency. It doesn't tag for ages:

https://github.com/ederc/GroebnerBasis.jl/issues/41

has issues with updating compats on time, and is GPL. That should not be deep in the dependency tree, so at most support via an addon package or Requires.

[blegat]

Yes, the current approach is for other packages to have SemialgebraicSets/MultivariatePolynomials in their dependency, e.g., HomotopyContinuation has SemialgebraicSets it its dependency to implement its interface and DynamicPolynomials/TypedPolynomials have MultivariatePolynomials in their dependency. However SemialgebraicSets and MultivariatePolynomials have a lightweight dependency tree.