Multi-Dimensional Bisection Method (MDBM) 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 in cases, where the number of unknowns surpasses the number of equations.
This type of problems can be found in many different field of science, just to mention a few:
- differential geometry (isolines, isosurfaces in higher dimensions)
- analysis of linkages (mechanical: workspace of robots)
- stability computation, stabilyzability diagrams
This method is an alternative to the contour plot or to isosurfaces in higher dimension, however, it has as the advantage of being able to handle multiple functions at once.
In addition, it uses far less function evaluation than the brute-force approach, making it much faster and more memory efficient, especially for complex tasks.
The bisection method - or the so-called interval halving method - is one of the simplest root-finding algorithms which is used to find zeros of continuous non-linear functions. This method is very robust and it always tends to the solution if the signs of the function values are different at the borders of the chosen initial interval.
Geometrically, root-finding algorithms of f(x)=0 find one intersection point of the graph of the function with the axis of the independent variable. In many applications, this 1-dimensional intersection problem must be extended to higher dimensions, e.g.: intersections of surfaces in a 3D space (volume), which can be described as a system on non-linear implicit equations. In higher dimensions, the existence of multiple solutions becomes very important, since the intersection of two surfaces can create multiple intersection curves.
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).
The software in this ecosystem was developed as part of academic research. If you use the MDBM.jl package as part of your research, teaching, or other work, I would be grateful if you could cite my corresponding publication: https://pp.bme.hu/me/article/view/1236/640
using MDBM using PyPlot; pygui(true);
Computation of a circle section defined by the implicit equation
foo(x,y) == 0 and by the constraint
c(x,y) > 0
function foo(x,y) x^2.0+y^2.0-2.0^2.0 end function c(x,y) #only the c>0 domain is analysed x-y end ax1=Axis([-5,-2.5,0,2.5,5],"x") # initial grid in x direction ax2=Axis(-5:2:5.0,"y") # initial grid in y direction mymdbm=MDBM_Problem(foo,[ax1,ax2],constraint=c) iteration=5 #number of refinements (resolution doubling) solve!(mymdbm,iteration) #points where the function foo was evaluated x_eval,y_eval=getevaluatedpoints(mymdbm) #interpolated points of the solution (approximately where foo(x,y) == 0 and c(x,y)>0) x_sol,y_sol=getinterpolatedsolution(mymdbm) fig = figure(1);clf() scatter(x_eval,y_eval,s=5) scatter(x_sol,y_sol,s=5)
Perform the line connection
myDT1=connect(mymdbm); for i in 1:length(myDT1) dt=myDT1[i] P1=getinterpolatedsolution(mymdbm.ncubes[dt],mymdbm) P2=getinterpolatedsolution(mymdbm.ncubes[dt],mymdbm) plot([P1,P2],[P1,P2], color="k") end
Intersection of two sphere in 3D
using LinearAlgebra axes=[-2:2,-2:2,-2:2] fig = figure(3);clf() #Sphere1 fS1(x...) = norm(x,2.0)-1 Sphere1mdbm=MDBM_Problem(fS1,axes) solve!(Sphere1mdbm,4) a_sol,b_sol,c_sol=getinterpolatedsolution(Sphere1mdbm) plot3D(a_sol,b_sol,c_sol,linestyle="", marker=".", markersize=1); #Sphere2 fS2(x...) = norm(x .- 0.5, 2.0) -1.0 Sphere2mdbm=MDBM_Problem(fS2,axes) solve!(Sphere2mdbm,4) a_sol,b_sol,c_sol=getinterpolatedsolution(Sphere2mdbm) plot3D(a_sol,b_sol,c_sol,linestyle="", marker=".", markersize=1); #Intersection fS12(x...) = (fS1(x...), fS2(x...)) Intersectmdbm=MDBM_Problem(fS12,axes) solve!(Intersectmdbm,6) a_sol,b_sol,c_sol=getinterpolatedsolution(Intersectmdbm) plot3D(a_sol,b_sol,c_sol,color="k",linestyle="", marker=".", markersize=2);
Mandelbrot set (example for a non-smooth problem)
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)
I am an assistant professor at the Budapest University of Technology and Economics, at the Faculty of Mechanical Engineering, the Department of Applied Mechanics.
During my studies and research, I have to determine stability charts of models described by delayed differential equations, which are typically formed as a "3 parameter / 2 implicit equation" problem. I have faced the difficulty that there is no applicable solution in any available software (e.g.: Mathematica, Matlab,...) which could easily be used in engineering problems.
Due to this reason, I have started to develop the Multi-Dimensional Bisection Method since 2006 in Matlab, and I have been improving it since, by adding new features from time to time.
I hope that others also find this package a useful tool in their work.
Best regards, Dr. Daniel Bachrathy
MDBM is a fiscally sponsored project of NumFOCUS, a nonprofit dedicated to supporting the open source scientific computing community.