How to solve two symbolic nonlinear equations in Julia?

Specific Domains Optimization (Mathematical)
1

Sorry if that seems an easy question but I searched for a simple answer in many packages without success. I have two equations like this:

𝐹_1: 𝑥^2+𝑥𝑦+\sin⁡(𝑦)=3
𝐹_2: 𝑥^3−4𝑥+\sin⁡(2𝑦)=0

and I need to solve them symbolically without giving initial guess or providing Jacobian or anything.

In MATLAB, for example, the following two lines do the job:

syms x y;
[sol_x,sol_y] = solve(x^2+x*y+sin(y) == 3, x^3-4*x+sin(2*y) == 0)

to get the answer:

sol_x = -2.0959358438825634582176494582886
sol_y = 1.0869388446654145864829239365948

Can I do something similar in Julia? I already searched in SymEngine.jl, Symata.jl, SymPy.jl, LinearExpressions.jl, etc but didn’t find a suitable answer.

1 Like
2

This kind of system of equations cannot be solved symbolically, you need a numerical approximation for this example.

The solution to that system is transcendental and not algebraic, so you can’t algebraically solve it.

Also, another symbolic package is Reduce.jl

4 Likes
3

Yeah, Matlab actually prints:

Warning: Cannot solve symbolically. Returning a numeric approximation instead. 
> In solve (line 303)

For numerical solution in Julia, you could use JuMP:

using JuMP, Ipopt
m = Model(solver=IpoptSolver())
@variable(m, x)
@variable(m, y)
@NLconstraint(m, x^2+x*y+sin(y) == 3)
@NLconstraint(m, x^3-4*x+sin(2*y) == 0)
solve(m)

after which you can get the solution found by Ipopt using

julia> getvalue.((x, y))
(2.0914461063943035, -0.4493417648042344)

Note that your set of equations has multiple solutions.

4 Likes
4

Or just use a tool made for this like NLsolve.jl. That would use a standard trust-region dogleg, Newton method, etc. based on your choice. I’m not sure exactly what IPOPT would be doing in this case. Is it equivalent to Newton or something?

1 Like
5

True, NLSolve.jl is the tool specifically for nonlinear equations root-finding. With JuMP, the root-finding is reformulated as an optimization problem with non linear constraints and IPOPT is used.

6

Or (overkill) use interval methods:

julia> using IntervalRootFinding, IntervalArithmetic, StaticArrays

julia> function f(z)
       x, y = z
       SVector(x^2+x*y+sin(y)-3, x^3-4*x+sin(2*y))
       end
f (generic function with 1 method)

julia> rts = roots(f, IntervalBox(-10..10, 2))
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([2.09144, 2.09145] × [-0.449342, -0.449341], :unique)
 Root([-2.09594, -2.09593] × [1.08693, 1.08694], :unique)
12 Likes
7

Just curious: Why is it overkill?

2 Likes
8

Now that I think about it, it probably isn’t. I am still adjusting to the power and simplicity of interval methods.

4 Likes
9

Interval methods are the only ones that can guarantee to have found all the roots in a given box, so in that sense they are not overkill.

4 Likes
10

I doubt you will be able to find symbolic solutions easily. Often such equations do not have closed-form solutions. Sometimes they can be solved with inverse functions, transformations or special functions as it is the case with the Lambert W function.

11

I thought it would be fun to check the solution graphically… using Wolfram Alpha (does Julia have an implicit plot function?)

image
image684×665 36 KB

2 Likes
12

You can make such plots within IntervalConstraintProgramming or ImplicitEquations:

using ImplicitEquations
using Plots
f(x,y) = x^2 + x*y + sin(y)
g(x,y) = x^3 - 4x + sin(2y)

plot(Eq(f,3), xlim=(-3pi,3pi), ylim=(-10,10))
plot!(Eq(g, 0), xlim=(-3pi,3pi), ylim=(-10,10))
5 Likes
13

Thanks! I’ll check it out!

14

Hm… I added package ImplicitEquations (Pkg.add("ImplicitEquations")), and tried to update the system (Pkg.update()), and got loads of error messages…

julia> Pkg.update()
INFO: Updating METADATA...
ERROR: METADATA cannot be updated. Resolve problems manually in C:\Users\username\.julia\v0.6\METADATA.
GitError(Code:EMERGECONFLICT, Class:Checkout, 425 conflicts prevent checkout)
macro expansion at .\libgit2\error.jl:99 [inlined]...

Anyway, I’m using Julia 0.6.3 on Win10 (64bit).

Another thing… if/when I get it to work…

  • If the implicit plot leads to a closed region, e.g., as in stability regions for Runge Kutta methods, is it possible to fill the interior or exterior of such a region by some filling color?
15

Hopefully it works out, if it does the answer is there is no filling support (though you can express inequalities in your functions, so you may be able to get what you want). The IntervalConstraintProgramming package might,you should check that out if it is a key feature.

16

The error means that for some reason you have local changes in METADATA, which stores information about package versions. METADATA is a git repository that Julia has cloned for you at C:\Users\username\.julia\v0.6\METADATA (from https://github.com/JuliaLang/METADATA.jl). cd into it, and then try git status to see what’s going on, probably git merge --abort to abort the merge, after which you can get rid of any local changes using git stash (undoable using git stash apply if you need it).

1 Like
17

Note that if you extend the domains of x and y, you will find more and more solutions (as I first observed using IntervalRootFinding.jl).

Here is the plot using IntervalConstraintProgramming. It shows how to plot filled regions, and how to easily intersect regions.

using IntervalArithmetic, IntervalConstraintProgramming, StaticArrays
using Plots

C1 = @constraint x^2 + x*y + sin(y) - 3 <= 0
C2 = @constraint x^3 - 4x + sin(2y) <= 0

X = IntervalBox(-20..20, 2)

paving1 = pave(C1, X, 0.1)
paving2 = pave(C2, X, 0.1)

plot( paving1.inner, lw=0, leg=false, aspect_ratio=1)
plot!(paving2.inner, lw=0)

plot!(paving.boundary, lw=1, lc=:blue, alpha=1, label="")
plot!(paving2.boundary, lw=1, lc=:red, alpha=1, label="")
#plot!(paving2.inner)

C3 = @constraint x^2 + x*y + sin(y) - 3 == 0
C4 = @constraint x^3 - 4x + sin(2y) == 0

C5 = C3 ∩ C4

paving3 = pave(C5, X, 0.1)

plot!(paving3.boundary, lc=:black, lw=5)

The resulting plot looks like this:

50%20AM
50%20AM944×882 116 KB

12 Likes
18

By playing with the tolerance, you can find all roots with the following code (using @Tamas_Papp’s definition of f). I would say that interval methods are not overkill – which other method would be able to find all of these roots, and guarantee that they are indeed all of the roots in that box?

julia> rts = roots(f, IntervalBox(-100..100, 2), 1e-14)
114-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([0.0404036, 0.0404037] × [98.879, 98.8791], :unique)
 Root([0.03142, 0.0314201] × [97.4523, 97.4524], :unique)
 Root([0.0208864, 0.0208865] × [95.7767, 95.7768], :unique)
 Root([0.0311381, 0.0311382] × [94.3102, 94.3103], :unique)
 Root([0.0431404, 0.0431405] × [92.5903, 92.5904], :unique)
 Root([0.033631, 0.0336311] × [91.1736, 91.1737], :unique)
 Root([0.0223543, 0.0223544] × [89.4906, 89.4907], :unique)
 Root([0.0333079, 0.033308] × [88.0313, 88.0314], :unique)
 Root([0.0462746, 0.0462747] × [86.3007, 86.3008], :unique)
 Root([0.0361763, 0.0361764] × [84.8955, 84.8956], :unique)
 ⋮
 Root([-0.0345102, -0.0345101] × [-84.8923, -84.8922], :unique)
 Root([-0.0231685, -0.0231684] × [-86.3475, -86.3474], :unique)
 Root([-0.0348574, -0.0348573] × [-88.0346, -88.0345], :unique)
 Root([-0.0446527, -0.0446526] × [-89.4457, -89.4456], :unique)
 Root([-0.0321866, -0.0321865] × [-91.1708, -91.1707], :unique)
 Root([-0.0215955, -0.0215954] × [-92.6338, -92.6337], :unique)
 Root([-0.0324881, -0.032488] × [-94.313, -94.3129], :unique)
 Root([-0.0417273, -0.0417272] × [-95.7348, -95.7347], :unique)
 Root([-0.0301559, -0.0301558] × [-97.4499, -97.4498], :unique)
 Root([-0.0202226, -0.0202225] × [-98.9197, -98.9196], :unique)

julia> all([rt.status for rt in rts] .== :unique)
true

By the way, for simpler functions you can even guarantee that you have found all (real) roots.
Using the master branch of IntervalRootFinding.jl, we have the following using an infinite initial box in 2D.
This is certainly impossible with any method that is not set-based.

julia> f(x, y) = SVector(x^2 + y^2 - 1, y - x)
f (generic function with 1 method)

julia> f(X) = f(X...)
f (generic function with 2 methods)

julia> rts = roots(f, IntervalBox(-∞..∞, 2))
2-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([0.707106, 0.707107] × [0.707106, 0.707107], :unique)
 Root([-0.707107, -0.707106] × [-0.707107, -0.707106], :unique)
9 Likes
19

I should emphasise that this plot is guaranteed, in the sense that what we are seeing are sets that are guaranteed to contain the true result.

1 Like
20

This is so cool…!

2 Likes