(feat) add falsi and bisection methods

Author: kanav99

Project.toml

@@ -0,0 +1,4 @@
+name = "NonlinearSolve"
+uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+authors = ["Kanav Gupta <[email protected]>"]
+version = "0.1.0"

src/NonlinearSolve.jl

@@ -0,0 +1,8 @@
+module NonlinearSolve
+
+  include("utils.jl")
+  include("bisection.jl")
+  include("falsi.jl")
+
+  export bisection, falsi
+end # module

src/bisection.jl

@@ -0,0 +1,71 @@
+"""
+  bisection(f, tup ; maxiters=1000)
+
+Uses bisection method to find the root of the function `f` between a tuple `tup` of values.
+"""
+function bisection(f, tup ; maxiters=1000)
+  x0, x1 = tup
+  fx0, fx1 = f(x0), f(x1)
+  fx0x1 = fx0 * fx1
+  fzero = zero(fx0x1)
+
+  (fx0x1 > fzero) && error("Non bracketing interval passed in bisection method.")
+  # NOTE: fx0x1 = 0 can mean that both fx0 and fx1 are very small and multiplication of them 
+  # could be less than the smallest float, hence could be zero.
+
+  fx0 == fzero && return x0 # should replace with some tolerance compare i.e. ≈
+  fx1 == fzero && return x1
+
+  left = x0
+  right = x1
+
+  iter = 0
+  while true
+    iter += 1
+    
+    if iter == maxiters
+      return left
+    end
+
+    fl = f(left)
+    fr = f(right)
+
+    fl * fr >= fzero && error("Bracket became non-containing in between iterations. This could mean that"
+    + " input function crosses the x axis multiple times. Bisection is not the right method to solve this.")
+
+    mid = (left + right) / 2.0
+    fm = f(mid)
+    if iszero(fm)
+      # we are in the region of zero, inner loop
+      right = mid
+      while true
+        iter += 1
+        
+        if iter == maxiters
+          return left
+        end
+        
+        mid = (left + right) / 2.0
+        (left == mid || right == mid) && return left
+        fm = f(mid)
+
+        if iszero(fm)
+          if !iszero(f(prevfloat_tdir(mid, x0, x1)))
+            return prevfloat_tdir(mid, x0, x1)
+          end
+          right = mid
+        else
+          left = mid
+        end
+
+      end
+    end
+
+    (left == mid || right == mid) && return left
+    if sign(fm) == sign(fl)
+      left = mid
+    else
+      right = mid
+    end
+  end
+end

src/falsi.jl

@@ -0,0 +1,74 @@
+"""
+  falsi(f, tup ; maxiters=1000)
+
+Uses Regula-Falsi method to find the root of the function `f` between a tuple `tup = (x0, x1)` of values.
+It doesn't find the exact value at which f(x) = 0, but finds an x where the next float in the direction of `x1`
+gives the evaluation of `f` to be zero. If `f` is zero for a region of x, it returns the left-most (in the 
+direction of x0) such value.
+"""
+function falsi(f, tup ; maxiters=1000)
+  x0, x1 = tup
+  fx0, fx1 = f(x0), f(x1)
+  fx0x1 = fx0 * fx1
+  fzero = zero(fx0x1)
+
+  (fx0x1 > fzero) && error("Non bracketing interval passed in bisection method.")
+  # NOTE: fx0x1 = 0 can mean that both fx0 and fx1 are very small and multiplication of them 
+  # could be less than the smallest float, hence could be zero.
+
+  fx0 == fzero && return x0 # should replace with some tolerance compare i.e. ≈
+  fx1 == fzero && return x1
+
+  left = x0
+  right = x1
+
+  iter = 0
+  while true
+    iter += 1
+    
+    if iter == maxiters
+      return left
+    end
+
+    fl = f(left)
+    fr = f(right)
+
+    fl * fr >= fzero && error("Bracket became non-containing in between iterations. This could mean that"
+    + " input function crosses the x axis multiple times. Bisection is not the right method to solve this.")
+
+    mid = (fr * left - fl * right) / (fr - fl)
+    fm = f(mid)
+    if iszero(fm)
+      # we are in the region of zero, inner loop
+      right = mid
+      while true
+        iter += 1
+
+        if iter == maxiters
+          return left
+        end
+        
+        mid = (left + right) / 2.0
+        (left == mid || right == mid) && return left
+        fm = f(mid)
+
+        if iszero(fm)
+          if !iszero(f(prevfloat_tdir(mid, x0, x1)))
+            return prevfloat_tdir(mid, x0, x1)
+          end
+          right = mid
+        else
+          left = mid
+        end
+
+      end
+    end
+
+    (left == mid || right == mid) && return left
+    if sign(fm) == sign(fl)
+      left = mid
+    else
+      right = mid
+    end
+  end
+end

src/utils.jl

@@ -0,0 +1,9 @@
+"""
+  prevfloat_tdir(x, x0, x1)
+
+Move `x` one floating point towards x0.
+"""
+function prevfloat_tdir(x::T, x0::T, x1::T)::T where {T}
+  x1 > x0 : prevfloat(x) : nextfloat(x)
+end
+