SymbolicUtils.jl icon indicating copy to clipboard operation
SymbolicUtils.jl copied to clipboard

Published 20 hours ago • verified publisher icon JuliaSymbolics

Equation solving

Open ChrisRackauckas opened this issue 5 years ago • 9 comments

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

ChrisRackauckas avatar May 05 '20 21:05 ChrisRackauckas

  • [ ] 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

dlfivefifty avatar May 06 '20 20:05 dlfivefifty

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

ChrisRackauckas avatar May 06 '20 20:05 ChrisRackauckas

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]}

dlfivefifty avatar May 07 '20 08:05 dlfivefifty

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 avatar Jul 27 '21 14:07 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.

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

anandijain avatar Jul 27 '21 14:07 anandijain

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)}

YingboMa avatar Jul 27 '21 18:07 YingboMa

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.

blegat avatar Jul 28 '21 08:07 blegat

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.

ChrisRackauckas avatar Jul 28 '21 11:07 ChrisRackauckas

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.

blegat avatar Jul 28 '21 12:07 blegat