# Hand-crafted Root Solver

This week our team was working on a problem where we needed to find roots of some fairly complex single-variable equations. Initially, we looked to the Ruby Standard Library which has a method, `#nlsolve`, that is designed to find roots. This solver ran a bit slow and got hung up on some initial guesses (it did the fridge stare). Looking at the implementation, we saw that this method was generalized to solve multi-dimensional systems which entailed finding the Jacobian matrix of the system and handing it off to an LU solver. This was overkill for our needs. We knew quite a bit about our equations. We knew the general shape of each curve, and we knew they each had a single root. There were also some constraints on the curve that allowed us to make a decent initial guess. Solving this problem wouldn’t be too complicated, so we decided to write our own root solver for a single equation, streamlining it to our needs. Here’s the secret sauce of how we did it:

The Bisection method is like a binary search for your root. You start with a high and low guess for a root. You choose high and low guesses that evaluate to values with opposite signs, meaning somewhere in between the values the curve crosses the x-axis. The method uses the average of the high and low guesses as the new guess. If the new guess is negative, it becomes the new low, and if it is positive, it becomes the new high. The method continues cutting its searching scope in half until you get close to a root or it determines there is no root in the scope. This method works well to converge quickly, even if you don’t have a great initial guess. It looks a little something like this:

```
def solve(f = @f, low = @low, high = @high, tol = @tol, n = @n)
x = (high + low) / 2
y = f.call(x)

if y.abs < tol #successfully found root
x
elsif n <= 0   #break out of solver after so many iterations
x
else #narrow window of searching by half
if y > 0
high = x
else
low  = x
end
solve(f, low, high, tol, n - 1)
end
end
```

For our problem, we use Bisection as a way to find a better first guess to pass off to Newton.

Newton’s method is an iterative approach that starts with an initial guess for a root, `x0`, and refines the guess by traveling along the curve based on its slope, where the new guess is defined as `x = x0 - f(x0)/f'(x0)`. This method is good for finding a root when you have a good initial guess. However, it can become unstable if it gets close to a horizontal slope (which is like dividing by zero in the equation above). There’s also the rare possibility that you could get stuck in an endless loop, where one guess cycles you back to a previous guess. Here’s a look at the core of our Newton’s method.

```
def solve(f = @f, x0 = @x0, tol = @tol, n = @n, eps = @eps)
y0 = f.call(x0)
y_prime = (f.call(x0 + eps) - y0) / eps

raise NonconvergenceError.new if y_prime.abs < tol

x = x0 - y0 / y_prime
y = f.call(x)

if y.abs < tol #successfully found root
x
elsif n <= 0   #not within threshold after n iterations.
x
else
solve(f, x, tol, n - 1, eps)
end
end
```

We packaged up our root solving methods in the aptly named gem, root_solver. The gem offers three methods for finding roots: Bisection, Newton, and a combination of the two that starts with Bisection then hands it off to Newton. Each requires a bit of configuration, which you can customize for a specific equation. The best thing you can do for this root solver is to choose your initial guess well. In addition, you can adjust the tolerance for root equality and the maximum number of iterations to run. This solver will return a value after a configured number of iterations, so the client may need to verify the solution if the number of iterations is low.

Serve It Up
All three root solvers in the gem respond to `#solve`, so that a client won’t have to know which root solving method it is using. The root solver object expects a callable (proc, lambda, PORO responding to call). The callable should take one argument, our `x`, and the callable should evaluate our function at `x`. Example callable:

```f = Proc.new { |x| x ** 2 - 10 } # Function for which we wish to find roots
l = 0                            # An educated guess for the lower limit
h = 5                            # An educated guess for the upper limit
t = 1e-3                         # How accurate do we want to be

RootSolver::BisectionNewton.new(f, l, h, t).solve
=> 3.1622776375625614

```

Share and Enjoy
We were able to speed up our calculations by an order of magnitude, just by simplifying the root solver to the special, but most common case, of a single equation.  The method in the Ruby Standard Library generalizes for a multi-dimensional system of equations, but it is likely that you just need a good, old-fashioned, single dependent variable Newton or Bisection solver like we provide in this gem.  If that’s your appetite, please, help yourself to our  root_solver, and enjoy.