# How to solve two symbolic nonlinear equations in Julia?

**URL:** <https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916>\
**Category:** Optimization (Mathematical)\
**Created:** [June 24, 2018, 10:58pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916 "2018-06-24T22:58:51Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [June 24, 2018, 10:58pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/1 "2018-06-24T22:58:51Z")

</div>

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:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [June 24, 2018, 11:44pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/2 "2018-06-24T23:44:30Z")

</div>

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

> **[GitHub - chakravala/Reduce.jl: Symbolic parser for Julia language term...](https://github.com/chakravala/Reduce.jl)**
>
> Symbolic parser for Julia language term rewriting using REDUCE algebra - GitHub - chakravala/Reduce.jl: Symbolic parser for Julia language term rewriting using REDUCE algebra

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [June 25, 2018, 5:49am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/3 "2018-06-25T05:49:33Z")

</div>

Yeah, Matlab actually prints:

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

```

For numerical solution in Julia, you could use JuMP:

```julia
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
julia> getvalue.((x, y))
(2.0914461063943035, -0.4493417648042344)

```

Note that your set of equations has multiple solutions.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 25, 2018, 6:30am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/4 "2018-06-25T06:30:47Z")

</div>

> [@tkoolen](#):
>
> For numerical solution in Julia, you could use JuMP:

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?

---

<div class="post-metadata">

**Author:** ![HCrepaz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hcrepaz/32/4531_2.png) [@HCrepaz](https://discourse.julialang.org/u/HCrepaz)\
**Post date:** [June 25, 2018, 6:47am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/5 "2018-06-25T06:47:26Z")

</div>

> [@ChrisRackauckas](#):
>
> Or just use a tool made for this like NLsolve.jl

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 25, 2018, 6:47am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/6 "2018-06-25T06:47:32Z")

</div>

Or (overkill) use interval methods:

```julia
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)  

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [June 25, 2018, 6:58am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/7 "2018-06-25T06:58:20Z")

</div>

Just curious: Why is it overkill?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 25, 2018, 7:02am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/8 "2018-06-25T07:02:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 25, 2018, 11:59am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/9 "2018-06-25T11:59:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![HCrepaz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hcrepaz/32/4531_2.png) [@HCrepaz](https://discourse.julialang.org/u/HCrepaz)\
**Post date:** [June 25, 2018, 1:04pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/10 "2018-06-25T13:04:03Z")

</div>

> [@Seif\_Shebl](#):
>
> and I need to solve them symbolically without giving initial guess or providing Jacobian or anything.

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.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [June 25, 2018, 5:56pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/11 "2018-06-25T17:56:41Z")

</div>

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/d/cd8ce20948c26eae0ce6c67f73fefb107a7c9a40.png)

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 25, 2018, 6:56pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/12 "2018-06-25T18:56:48Z")

</div>

You can make such plots within `IntervalConstraintProgramming` or `ImplicitEquations`:

```julia
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))

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [June 25, 2018, 8:26pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/13 "2018-06-25T20:26:01Z")

</div>

Thanks! I’ll check it out!

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [June 25, 2018, 8:49pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/14 "2018-06-25T20:49:43Z")

</div>

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [June 25, 2018, 8:55pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/15 "2018-06-25T20:55:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [June 26, 2018, 7:57am UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/16 "2018-06-26T07:57:56Z")

</div>

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](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).

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 26, 2018, 2:32pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/17 "2018-06-26T14:32:22Z")

</div>

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.

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/7/e/7e33261e39ba99adf8b8f7910709ca3dfd021094.jpeg)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 26, 2018, 2:39pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/18 "2018-06-26T14:39:13Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 26, 2018, 2:39pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/19 "2018-06-26T14:39:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [June 26, 2018, 2:51pm UTC](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916/20 "2018-06-26T14:51:39Z")

</div>

This is so cool…!

[Next page](https://discourse.julialang.org/t/how-to-solve-two-symbolic-nonlinear-equations-in-julia/11916.md?page=2)