# \[ANN\] MDBM.jl - A root-finding package for system of non-linear equations

**URL:** https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141
**Category:** Package Announcements
**Tags:** roots
**Created:** [March 21, 2019, 1:14pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141 "2019-03-21T13:14:39Z")
**Posts on this page:** 13
**Page:** 1

<div class="post-metadata">

### Author: ![bachrathyd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bachrathyd/32/7545_2.png) [@bachrathyd](https://discourse.julialang.org/u/bachrathyd)
#### Post date: [March 21, 2019, 1:14pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/1 "2019-03-21T13:14:39Z")

</div>

Multi-Dimensional Bisection Method

[MDBM.jl](https://github.com/bachrathyd/MDBM.jl) is a package is an efficient and robust root-finding algorithm, which can be used to determine whole high-dimensional submanifolds (points, curves, surfaces…) of the roots of implicit non-linear equation systems, especially of the cases, where the number of unknowns surpasses the number of equations.

It uses far less function evaluation than the brute-force approach, making it much faster and more memory efficient, especially for complex tasks.

MDBM algorithm can handle:

- multiple solutions
- arbitrary number of variable (typically: 3-6)
- arbitrary number of implicit equations
- multiple constraints in the parameter space
- degenerated functions

while providing the gradients of the equations for the roots by means of first order interpolation (and convergence rate).

This method can be used to create contour plots or isosurfaces in higher dimension, furthermore, it has as the advantage of being able to handle multiple functions at once.

It is a Julia implementation of a [Matlab package](https://www.mathworks.com/matlabcentral/fileexchange/69414-multi-dimensional-bisection-method).  
Please have a look at [MDBM.jl](https://github.com/bachrathyd/MDBM.jl) and at the examples. Feel free to suggest improvements, design choices and submit PR.

I hope you will find the algorithms useful in your own work.

Best wishes,  
Daniel Bachrathy

---

<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: [March 21, 2019, 2:20pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/2 "2019-03-21T14:20:35Z")

</div>

Neat! Can you elaborate a bit about the requirements for the functions it can handle? Eg I understand that they need be continuous, but do they need to be ADable or evaluating reals to reals is sufficient?

Also, can your provide some intuition about how the solution time scales with the dimensionality of the unknowns and the number of equations?

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [March 21, 2019, 4:54pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/3 "2019-03-21T16:54:01Z")

</div>

Very cool! In addition to Tamas’ questions would you mind giving a bit of a run-down of how this compares to NLsolve.jl, ie. relative strengths and weaknesses?

---

<div class="post-metadata">

### Author: ![bachrathyd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bachrathyd/32/7545_2.png) [@bachrathyd](https://discourse.julialang.org/u/bachrathyd)
#### Post date: [March 21, 2019, 5:13pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/4 "2019-03-21T17:13:23Z")

</div>

There is no any restriction for the function, because actually the method finds the ‘sing-change’ and not the ‘root’.  
So in case of f=1/x, the solution provided by the method will converge to 0 (despite to the ill conditioned interpolation).

Even, degenerated functions can be used (see the example at the end)

So it is sufficient that f evaluates reals^n to reals^m (if m\<=n).

The complexity of the method is proportional to the (fractal) dimension of the object  
e.g:  
3 param & 2 equation → curves: 1D object, so the computation time doubles (2^1) every iteration steps  
4 param & 2 equation → surface: 2D object, so the computation time quadruple every iteration steps

I forgot to include, but I gave a talk at [JuliaCon2018/MDBM](https://www.youtube.com/watch?v=AKCZmWStRU0)

The computation needs are explained at [8:09](https://youtu.be/AKCZmWStRU0?t=489).

Degenerated example:

```julia
ax1=Axis(-3.0:3.0,"x")
ax2=Axis(-3.0:3.0,"b")

function foo(x,y)
    if norm([x,y])<1.0
        0.0
    else
        x-y
    end
end

mymdbm=MDBM_Problem(foo,[ax1,ax2])#ax3
iteration=5 #number of refinements (resolution doubling)
solve!(mymdbm,iteration,interpolationorder=1);

#evaluated points
x_eval,y_eval=getevaluatedpoints(mymdbm)

#solution points
x_sol,y_sol=getinterpolatedsolution(mymdbm)

fig = figure(1);clf()
scatter(x_eval,y_eval,s=2)
scatter(x_sol,y_sol,s=4);

```

![degen](https://global.discourse-cdn.com/julialang/original/3X/0/0/00ff32e61d8c399375803e09bd91a76c3e1f67e2.png)

---

<div class="post-metadata">

### Author: ![bachrathyd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bachrathyd/32/7545_2.png) [@bachrathyd](https://discourse.julialang.org/u/bachrathyd)
#### Post date: [March 21, 2019, 5:26pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/5 "2019-03-21T17:26:50Z")

</div>

I did not used the NLsolve.jl, however, it seems that is provides a single solution only.  
The great strengths of my method that is can solve problem with more parameter than equations  
e.g.:  
f(x,y)=x^2+y^2-1.0  
which is typical in many problem (computation of a chart, a bifurcation curve,…).  
Furthermore, it can handle non-smooth functions (see the previous comment).  
This type of problem, where the existence of multiple solution is important (especially as a function of one or two parameter) the MDBM is excellent.

Compared it to a continuation method, the MDBM is automatic (e.g.: no need to start a new branch), and it can find closed curves (isola) which is hard to detect by continuation.

The drawbacks of MDBM

- that is has only linear convergence (but it is not a big problem in chart plotting)
- it can also miss small closed curves (isola) if the initial resolution is not fine enough
- handling more the 6 parameter or computing higher order objects (surfaces, volumes…) can be very memory & time consuming

---

<div class="post-metadata">

### Author: ![bachrathyd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bachrathyd/32/7545_2.png) [@bachrathyd](https://discourse.julialang.org/u/bachrathyd)
#### Post date: [March 21, 2019, 5:28pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/6 "2019-03-21T17:28:37Z")

</div>

p.s.: we have found a small bug, which is fixed in the new 0.1.3 release!

---

<div class="post-metadata">

### Author: ![StevenSiew](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevensiew/32/218393_2.png) [@StevenSiew](https://discourse.julialang.org/u/StevenSiew)
#### Post date: [March 23, 2019, 6:06am UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/7 "2019-03-23T06:06:56Z")

</div>

Thank you for creating this module. It is much appreciated.

---

<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: [March 23, 2019, 7:37am UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/8 "2019-03-23T07:37:59Z")

</div>

Very nice! How does this compare to IntervalRootFinding.jl?

---

<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: [March 23, 2019, 7:38am UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/9 "2019-03-23T07:38:53Z")

</div>

By the way, I think it would be useful to change the package name to be more explicit and discoverable.

---

<div class="post-metadata">

### Author: ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)
#### Post date: [March 23, 2019, 8:18am UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/10 "2019-03-23T08:18:53Z")

</div>

Very cool! Look forward to trying it out on some problems I have!

---

<div class="post-metadata">

### Author: ![bachrathyd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bachrathyd/32/7545_2.png) [@bachrathyd](https://discourse.julialang.org/u/bachrathyd)
#### Post date: [March 25, 2019, 1:21pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/11 "2019-03-25T13:21:08Z")

</div>

I think that the `IntervalRootFinding.jl` is definitely better for simpler analytical expressions than the `MDBM.jl`, because it can provide guaranteed and unique solutions, while in the `MDBM` there is slight a chance to lose some solutions, however the usage of the `IntervalRootFinding.jl` is limited.

The `MDBM` is initially designed for stability chart computations which is based on the eigenvalues of a large matrix, for which the `IntervalRootFinding` would fail (see the first example below).  
So the `MDBM` can be used for broader problems, e.g.: eigenvalue based (1st example) and iteration based (2nd example) functions.  
Furthermore, it was used to ‘black-box’ functions where the output of a Finite Element Software (Abaqus) is used. The `MDBM` method was successfully used in automatically evaluated measurements to find the stability boundary of a controlled mechanical system, where noise was significant (see the third example with nois.

1# Egienvalue problem

```julia
using LinearAlgebra, Arpack
function fooeig(x)
       mymx=Matrix{Float64}(I,5,5);
       mymx[2,1]=x;
       mymx[1,2]=x+5;
       abs(eigs(mymx)[1][1])-3
end
eig_problem=MDBM_Problem(fooeig,[Axis(-10:10)])
solve!(eig_problem,10)
getinterpolatedsolution(eig_problem)

```

Roots:

```julia
1-element Array{Array{Float64,1},1}:
 [-5.70156, 0.701562]

```

#2. Iteration based problem:

```julia
function mandelbrot(x,y)    
    c=x+y*im
    z=Complex(zero(c))
    k=0
    maxiteration=1000
    while (k<maxiteration && abs(z)<4)
            z=z^2+c
            k=k+1
        end
    return abs(z)-2
end

Mandelbrotmdbm=MDBM_Problem(mandelbrot,[-5:2,-2:2])
solve!(Mandelbrotmdbm,9)
a_sol,b_sol=getinterpolatedsolution(Mandelbrotmdbm)
fig = figure(3);clf()
plot(a_sol,b_sol,linestyle="", marker=".", markersize=1)

```

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

#3. Function with noise

```julia
function foo(x,y)
    x^2.0+y^2.0-2.0^2.0+rand()
end

mymdbm=MDBM_Problem(foo,[-3.0:3.0,-3.0:3.0])
solve!(mymdbm,5)

x_sol,y_sol=getinterpolatedsolution(mymdbm)
fig = figure(1);clf()
scatter(x_sol,y_sol,s=4);

```

![noisy_circle](https://global.discourse-cdn.com/julialang/original/3X/4/8/4836b8e5d43b40358645f10340aecdd1e7ccd313.png)

---

<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: [March 25, 2019, 2:48pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/12 "2019-03-25T14:48:10Z")

</div>

Nice examples, thanks!

---

<div class="post-metadata">

### Author: ![cce](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cce/32/460_2.png) [@cce](https://discourse.julialang.org/u/cce)
#### Post date: [April 4, 2019, 1:29pm UTC](https://discourse.julialang.org/t/ann-mdbm-jl-a-root-finding-package-for-system-of-non-linear-equations/22141/13 "2019-04-04T13:29:37Z")

</div>

Dr. Bachrathy,

From the peanut gallery with pedantic thoughts… I comment on using acronym for the package: it could have been named `MultiDimensionalBisection.jl`. There’s no harm in a “using MultiDimensionalBisection” at the top of a file, and actually, it’s quite helpful/informative. On the other hand, `MDBM.jl` made me think it’s some sort of database or … well, acronyms are not as informative as longer, spelled-out package names. This makes it harder to search for on google as well. I don’t know if it’d even be possible to change the name, etc. Regardless, I love the examples: as someone who is not particularly numerically inclined they were quite illustrative.

Congragulations,  
Clark
