RootFinding GSOC project #1192

n0rbed · 2024-07-25T12:58:06Z

Hello! this is a PR concerning the GSOC project "RootFinding for Symbolics.jl" which you can check more details about here. The aim of this project is to add the solve function to Symbolics.jl for solving equations symbolically. The supervisors of this project are @sumiya11 and @shashi.

Feel free to try the solver out yourself. Comments and suggestions are very welcome.

There are 2 issues in this PR regarding licenses: (1) Nemo.jl may be not LGPL v3 compatible, we are working on it. (2) Groebner.jl is GPL, so it goes into an extension.

Prior to this PR, all work and progress was carried out here

Currently, there are 4 available solvers which are grouped under one function, solve.

solve_univar
solve_multipoly
solve_multivar
ia_solve

To give you an idea of what these solvers are capable of, here are some examples:

`solve_univar` (uses factoring and analytic solutions up to degree 4)

julia> expr = expand((x   3)*(x^2   2x   1)*(x   2))
6   17x   17(x^2)   7(x^3)   x^4

julia> solve(expr, x)
3-element Vector{Any}:
 -2//1
 -1//1
 -3//1

julia> solve(expr, x, true)
4-element Vector{Any}:
 -2//1
 -1//1
 -1//1
 -3//1

julia> f = expand((x^3 - 1) * (c*x^2   a*x   b) * (x^10 - 1))
b   a*x   c*(x^2) - b*(x^3) - a*(x^4) - c*(x^5) - b*(x^10) - a*(x^11) - c*(x^12)   b*(x^13)   a*(x^14)   c*(x^15)

julia> Symbolics.solve(f, x)
14-element Vector{Any}:
   (1//4)   (1//2)*(Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)) - 4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))) - Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))
   (1//4)   (1//2)*(-Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)) - 4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))) - Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))
   (1//4)   (1//2)*(Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))   Symbolics.ssqrt(-4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   -5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))   (1//64)*(-80   Symbolics.ssqrt(5120))))
   (1//4)   (1//2)*(Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))) - Symbolics.ssqrt(-4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   -5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))   (1//64)*(-80   Symbolics.ssqrt(5120))))
 -1
  1
   -(1//4)   (1//2)*(-Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))   Symbolics.ssqrt(-4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   -5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))   (1//64)*(-80   Symbolics.ssqrt(5120))))
   -(1//4)   (1//2)*(-Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))) - Symbolics.ssqrt(-4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   -5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)))))   (1//64)*(-80   Symbolics.ssqrt(5120))))
   -(1//4)   (1//2)*(Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)) - 4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))))   Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))
   -(1//4)   (1//2)*(-Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120)) - 4((5//16)   (1//128)*(-80   Symbolics.ssqrt(5120))   5 / (16Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))))   Symbolics.ssqrt((1//64)*(-80   Symbolics.ssqrt(5120))))
   (-a   Symbolics.ssqrt(a^2 - 4b*c)) / (2c)
   (-a - Symbolics.ssqrt(a^2 - 4b*c)) / (2c)
   (1//2)*(-1   Symbolics.ssqrt(-3))
   (1//2)*(-1 - Symbolics.ssqrt(-3))

# We can see here that it factors this poly to ((1//1)   x   x^2   x^3   x^4   x^5   x^6)(x-1//1). 
# One of these polys can be solved by analytic methods and the other is not, so the unsolvable
# factor is just returned in a roots_of struct.
julia> solve(x^7 - 1, x)
2-element Vector{Any}:
roots_of((1//1)   x   x^2   x^3   x^4   x^5   x^6)
 1//1

`solve_multivar` (uses Groebner basis and `solve_univar` to find roots)

julia> eqs = [x y^2 z, z*x*y, z 3x y]
3-element Vector{Num}:
 x   z   y^2
       x*y*z
  3x   y   z

julia> solve(eqs, [x,y,z])
3-element Vector{Any}:
 Dict{Num, Any}(z => 0//1, y => 0//1, x => 0//1)
 Dict{Num, Any}(z => 0//1, y => 1//3, x => -1//9)
 Dict{Num, Any}(z => -1//1, y => 1//1, x => 0//1)

`solve_multipoly` (uses GCD between the input polys)

julia> solve([x-1, x^3 - 1, x^2 - 1, (x-1)^20], x)
1-element Vector{Rational{BigInt}}:
 1

`ia_solve` (solving by isolation and attraction)

This function attempts to follow the PRESS system in [1]. It first checks the number of occurrences, then uses isolate or attract as needed.

attract(lhs, var) currently uses 4 techniques for attraction.

Log addition: log(f(x)) log(g(x)) => log(h(x))
Exponential simplification: a*b^(f(x)) c*d^(g(x)) => f(x) * log(b) - g(x) * log(d) log(-a/c).
Trig simplification: this bruteforces multiple trig identities and doesn't detect them before hand.
Polynomialization: as a last resort, attract attempts to polynomialize the expression. Say sin(x 2)^2 sin(x 2) 10 is converted to X^2 X 10.

julia> solve(2^(x 1)   5^(x 3), x)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 (-log(2)   3log(5) - log(complex(-1))) / (log(2) - log(5))

julia> solve(log(x 1) log(x-1), x)
2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 (1//2)*RootFinding.ssqrt(8.0)
 (-1//2)*RootFinding.ssqrt(8.0)

issue #866 😉

julia> solve(a*x^b   c, x)
((-c)^(1 / b)) / (a^(1 / b))

Issues solved

issue #1115

julia> Symbolics.solve(2x^4   7x^3   4x^2   7x   2, x)
4-element Vector{Any}:
       (1//4)*(-7   Symbolics.ssqrt(33))
       (1//4)*(-7 - Symbolics.ssqrt(33))
 0//1   1//1*im
 0//1 - 1//1*im

issue #961, similarly, issue #525 and #524 are also solved

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

julia> eqs = [x*y - 4, x   y - 4]
2-element Vector{Num}:
   -4   x*y
 -4   x   y

julia> Symbolics.solve(eqs, [x,y])
2-element Vector{Any}:
 Dict{Num, Any}(y => 2, x => 2)
 Dict{Num, Any}(y => 2, x => 2)

issue #424

julia> Symbolics.solve(b/a - 1, b)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
a

Difficulties faced

Whilst developing the polynomial solvers, complex numbers and domains were found to cause problems in multiple aspects of our results. First, sqrt needed another complex() inside of it if the integer is negative, say sqrt(-2) does not work, sqrt(complex(-2)) works. Cbrts however, dont support complex nums, so you'll have to do (complex(1))^(1/3) if you wanna compute the cuberoot of 1 0im. Substituting variables by complex numbers also caused bugs. To solve this, included in the file solve_helpers.jl and ia_helpers.jl are helper functions such as ssubs, ssqrt, scbrt, and slog which handle all these cases themselves so that the output is nice and tidy, and that the program at runtime does not need to evaluate whatever is in a square root for example to check if its inside our domain or not.

Documentation

There's an in depth documentation here which covers the first 3 solvers in depth since they were not mentioned here. Some parts of this documentation may be outdated though, so do be careful since functions may have been renamed or altered slightly.

References

[1]: R. W. Hamming, Coding and Information Theory, ScienceDirect, 1980.

codecov-commenter · 2024-07-25T13:03:35Z

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 89.63265% with 127 lines in your changes missing coverage. Please review.

Project coverage is 79.42%. Comparing base (79c4e92) to head (fe245e1).
Report is 493 commits behind head on master.

Files with missing lines	Patch %	Lines
src/solver/ia_main.jl	81.32%	31 Missing ⚠️
src/solver/polynomialization.jl	83.44%	25 Missing ⚠️
src/solver/solve_helpers.jl	79.16%	20 Missing ⚠️
src/solver/attract.jl	80.89%	17 Missing ⚠️
src/solver/main.jl	90.00%	12 Missing ⚠️
ext/SymbolicsGroebnerExt.jl	95.78%	7 Missing ⚠️
src/solver/nemo_stuff.jl	63.63%	4 Missing ⚠️
src/solver/univar.jl	95.65%	4 Missing ⚠️
src/solver/ia_helpers.jl	96.47%	3 Missing ⚠️
src/solver/preprocess.jl	98.50%	2 Missing ⚠️
... and 2 more

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1192       /-   ##
==========================================
  Coverage   77.07%   79.42%    2.34%     
==========================================
  Files          32       46       14     
  Lines        3533     4626     1093     
==========================================
  Hits         2723     3674      951     
- Misses        810      952      142

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Project.toml

src/Symbolics.jl

src/solver/RootFinding.jl

src/solver/main.jl

test/runtests.jl

sumiya11 · 2024-07-26T20:53:05Z

@n0rbed thank you very much for making the PR. I will add some comments.

ext/SymbolicsGroebnerExt.jl

src/solver/attract.jl

src/solver/coeffs.jl

src/solver/attract.jl

src/solver/ia_helpers.jl

src/solver/ia_main.jl

src/solver/main.jl

src/solver/solve_helpers.jl

src/solver/postprocess.jl

…loat as it is, also removed trig unused func

ext/SymbolicsGroebnerExt.jl

src/solver/polynomialization.jl

src/solver/solve_helpers.jl

…into RootFinding

RootFinding GSOC project

0d6395e

n0rbed added 2 commits July 25, 2024 16:03

spelling mistakes

3535116

detect difficult poly case

e5bee71

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

Project.toml Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

src/Symbolics.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

src/solver/RootFinding.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

src/solver/main.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

src/solver/main.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jul 26, 2024

View reviewed changes

test/runtests.jl Outdated Show resolved Hide resolved

added g(x) sqrt(f(x)) c solution in attract

9f3dcc7