One dimensional root finding algorithms

• Download source files - 258 Kb
• Download demo project - 114 Kb

Introduction

This article talks about root finding in one dimension. It discusses about some algorithms than can solve the problem f(x) = 0, where f is a given function and x is a real variable.

Background

Many algorithms can be found on the Internet, with their source code. My job here was just to translate them in C#, sometimes modify them a bit, and essentially embed them into a convenient framework, a la, C++ STL.

About root finding

How can a computer solve an equation like f(x)=0? As you can imagine, root finding algorithms don't "solve" the equation. They "only" provide (in the best case) one approximated solution, using iterative methods. Starting with a given interval, that is assumed to contain the solution, the algorihtm reduces at least by 2 (using the Bisection method) the length of the interval at each iteration. It stops when the length of the remaining interval reaches a pre-defined value (called Accuracy).

Basically, this works fine. But it can happen that the Accuracy can not be reached. To prevent the algorithm from iterating infinitely, you should prefer to stop it when the number of iterations reaches a certain value (called Iterations).

Now, you're almost ready to use the code. But you should be interested with a special improvement. Suppose you don't know the function f you're working with. How can you determine the range the algorithm has to start with? In this case, you have to use the OutwardBracketing. Given a growing factor, the algorithm will iteratively grow the initial range until the function f changes its sign between the bounds of the range. Note that because f is supposed to be continuous, this ensures that a solution lies in the grown interval.

Using the code

Using the library is very easy. The library provides six methods, among the most usual ones. Each one is embedded in a class, which itself inherited from the abstract RootFinder class, like it is shown on the picture:

The most difficult of your job is to choose the good method. This will be discussed later. Suppose you want to solve the equation f(x) = 0, where f is defined like this:

private double f (double x) {
    return 3x*x-1;
}

The two solutions lie in the interval [-1;+1]. Suppose you want the positive root. Don't forget that the algorithms provide a unique solution. If you choose a too big starting interval, the algorithm may go wrong, or give you the solution you don't want. So, you have to start with an interval like [0;1].

private double FindRoot (double x) {

     // Create the root finder object, that contains the algorithm
     BisectionRootFinder finder =
            new BisectionRootFinder (new UnaryFunction(f));

     // Define the accuracy you want for the root
     finder.Accuracy = 1.0E-04;

     // Prevent overflow
     finder.Iterations = 30;

     // Compute without bracketing outward
     double root = finder.Solve(0.0,1.0,false);

     return root;
}

The root r of f on [0;1] is r = sqrt(1/3). The Accuracy ensures you that the approximated root r' provided by the algorithm will verify :

Abs(r'-r) < 10e-4 (Accuracy)

except if the algorithm reaches the maximum of iterations you passed it (30 here).

Which algorithm may you choose?

If you have an expression of the derivative function df of f, use the Newton-Raphson method. It's stable and converges quickly. If you don't have an expression df of the derivative function of f, the answer is not easy. It depends, in fact, on what is the most important for you. The Bisection method will always lead you to a good approximation, but not as quickly as the others. The three methods False Position, Secant, and Ridder converge more quickly, but can fail in some cases. The Brent method seems so to be the most efficient method.

Points of Interest

The code source of the algorithms can be found everywhere on the Internet. Building a framework for one dimension functions is not easy in C#. In C++, you can create functors and arithmetic operators on template classes. In C#, you cannot. Okay, handling pointers on functions is very easy with delegates, but that's all. So, how can we build an efficient framework for functions in C#? To be honest, I think it is impossible. That's the reason why I don't provide such a framework. Creating an abstract class Function and derive it as needed is not necessarily the solution. What I propose in the downloadable package is a static class, named <c>Functions that can help you with handling the functions: addition, substraction, and composition.

In the demo package, you will find a mathematics formula parser. The original idea and the code comes from Marcin Cuprjak's parser. The Chart control I used to plot the graph of the functions is from JChampion, and can be downloadable here. The rich panels were designed by Prasad Khandekar.