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
`` js
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;
}
`

$3

|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

$3

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

$3

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)$.




$3


- 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.

$3

` js
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

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

$3

` js
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

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:
` js
function myZero(f, a, b, max_iter, yTol) {
const gen = zero_generator(a, b);

let result = gen.next();
for (let i=0; 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.

!Picture valid region for IQI

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