chandrupatla-zero

This module implements Chandrupatla's algorithm to find a root of a function f(). This is a little known fast and efficient algorithm, and in my opinion outperforms Brent's zero().

The module provides 2 functions: zero(), and zero_generator(). The function zero_generator() facilitates solving asynchronous functions, and provide a lot of flexibility (see example below). For typical use, the simpler zero() function should be sufficient.

Usage

Chandrupatla's zero(), like Brent's zero() retuns a zero x of the function f in the given interval [x1, x2], to within a tolerance. The tolerance is the maximum of:

1. 2 times smallest positive number greater than 0
   2 * Number.MIN_VALUE
2. 2 times the smallest step to the next double precision number from the maximum of the absolute values of the intervals edges
   2 ** (frexp(Math.min(Math.abs(x1), Math.abs(x2)))[1] + machep).
   The function frexp() is imported from the package locutus/c/math.
3. a user defined function options.userTolerance(x1, x2, f1, f2)
4. a calculation based on the optional parameters options.eps and options.dlt

To create a function userTolerance() similar to this calculation, used in the code

function userTolerance(x1, f1, x2, f2) {
    const xm = (Math.abs(f1) < Math.abs(f2) ? x1 : x2);
    return 2 * eps * Math.abs(xm) + 0.5 * dlt;
}

Parameters

Parameter	Description
f	the function f() for which the root is sought. Only for the utility function zero().
x1	the finite lowerbound of the interval in which to search the root
x2	the finite upperbound of the interval in which to search the root
options	the dictionary containing the optional parameters:
	- f1, f2: posibly already known function values for f(x1), and f(x2)
	- x3: the first X value to try, instead of the bisection value
	- userTolerance: a function taking 2 parameter to calculate an alternative tolerance
	- eps: a relative tolerance
	- dlt: an absolute tolerance

Functions

• zero_generator(a, b, opts), and
• zero(f, a, b, opts).

Yields and returns

The generator functions zero_generator() yields 'x' values to request function evaluations for those values. When the zero() is called, or the zero_generator() is 'done', they return an array with the following contents:

• res[0]: the x-value with the smallest absolute function value $f(x)$,
• res[1]: the corresponding function value $f(x)$,
• res[2]: the other x-value for the interval, and
• res[3]: its corresponding function value $f(x)$.

Notes

• f(x1) and f(x2) should have different signs.
• the algorithm ends when f(x) == 0 is found, or the interval [a, b] has been sufficiently shrunk to be sure the x is within the 2 times the tolerance. This says nothing about the value of f(x). For example the function x => x < 3 ? -1 : 2, will find a "root" close to 3, but f(x) will be -1.

Some Examples

The following examples shows simple usage. Given a function f(x) = x² - 5.

1. Simplest zero() example

const zeros = require('chandrupatla-zero');

console.log( zero(x => x * x - 5, 1, 4) );
// expect:
// [
//   2.23606797749979,
//   8.881784197001252e-16,
//   2.236067977499789,
//   -3.552713678800501e-15
]

with known f1=-4 and f2=11, and eps=Number.EPSILON

console.log( zero(x => x * x - 5, 1, 4, {
   f1: -4,
   f2: 11,
   eps: Number.EPSILON,
})[0].toFixed(9) );
// expect:
// 2.236067977

2. Simple zero_generator() example

const gen = zero_generator(-3, -2);

let result = gen.next();
while ( ! result.done ) {
    const x = result.value;
    const y = f(x);

    result = gen.next(y)    
}
console.log( result.value[0].toFixed(6), result.value[1] );
// expect:
// -2.236068 8.881784197001252e-16

3. Modified zero_generator() example

Often implementations have additional parameters for maximum number of iterations (max_iter), or an allowed tolerance for the root's function value (yTol). Chandrupatla's original code does not have these options, and I am fairly sure he would not agree with them. However, using the generator function, they are easily implemented, as the following example function shows:

function myZero(f, a, b, max_iter, yTol) {
    const gen = zero_generator(a, b);

    let result = gen.next();
    for (let i=0; i<max_iter && ! result.done; ++i) {
        const x = result.value;
        const y = f(x);
        if (Math.abs(y) < yTol) break;

        result = gen.next(y)    
    }

    return result.value;
}

Chandrupatla's Algorithm

The algorithm always works with 3 points. The points are numbered:

1. The point x1 is the middle point. At the start of the algorithm it is determined using bisection.
2. The point x2 is on the opposite side of the X-axis from the point x1. This is f(x_1) f(x_2) < 0.
3. The point x3 is the discarded point (on the same side of the X-axis as the point x1).

In each iteration of the algorithm we start by scaling points using
scale(x, y) = ((x - x2)/(x3-x2), (f(x)-f(x2))/(f(x3)-f(x2)))
This will scale (x2, f(x2)) to (0, 0), and scale (x3, f(x3)) to (1, 1).

Then we determine of the scaled point for (x1, f(x1)) is in the "valid region". The "valid region" is the region were the points combined with (0, 0) and (1, 1) will result in an Inverted Quadratic Interpolation (IQI) parabole whose top has a Y-value outside the interval [0, 1].

If the scaled point (x1, f(x1)) is inside the "valid region" we use IQI to determine the next points X-value. Otherwise the next points X-value is determined by bisection.

Further Documentation

The ultimate documentation on this algorithm can be found on Emeritus Professor T.R. Chandrupatla's article:

Chandrupatla, Tirupathi R., A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives. Advances in Engineering Software, 1997, 28 (3): 145–149. doi:10.1016/S0965-9978(96)00051-8. .

License

ISC