0

With root-finding algorithms, you're usually just given the x value of the position of a root. I want to get the two closest points on either side of the root where these two points would have opposite signs considering that's where the roots exist (where the y value goes from + to - or the other way around)

So given an x value that is the approximate result of a root within an interval, what's the best efficient method to find the nearest two points to that root that have opposite y values.

Currently, my approach is to keep widening the offset from the roots position on either side until the y values of these two points are opposite.

Also, maybe using the precision (like 1e-6) that was used to find the root could be helpful to determine how much of an offset is needed to get passed the actual root where the signs change!

enter image description here

This image is sort of what I want. Blue are the two points I want that are at either side of the approximate root (which is purple). Of course the closer they are the better (without using further more iterations to get a more precise root, but instead only use the root given - we can assume that it's close enough).

I know some algorithms like bisection give the last closest brackets that could be used as these points unless the 1st iteration found the root exactly at the midpoint and the brackets are still far apart. While I'm using ridders (regula falsi) method which has a bias so one bracket usually doesn't move.

Code

let x = 3.135, y = fn(x); /* x = approximate result of root, actual root is pi */
let xmin = 3.1, xmax = 3.2;
let offset = (xmax - xmin) * .0000000002;

function getOffsetPoints(x, y, offset, fn) {

        /* if root exactly 0 */
        /* offset either side by 1e-15 = .000000000000001! */
        if (y === 0) {
            return [
                [x - 1e-15, fn(x - 1e-15)],
                [x + 1e-15, fn(x + 1e-15)]
            ]
        }

        let dxLeft = x - offset,
            dyLeft = fn(dxLeft),
            dxRight = x + offset,
            dyRight = fn(dxRight);

        /* maybe 3 attempts is enough - don't want to over do it with computations */
        let i = 0;
        while (i <= 3 && Math.sign(dyLeft) * Math.sign(dyRight) !== -1) {
            offset *= 100;
            dxLeft = x - offset;
            dyLeft = fn(dxLeft);
            dxRight = x + offset;
            dyRight = fn(dxRight)
            i++;
        }

        /* return two points at either side of root */
        return [
            [dxLeft, dyLeft],
            [dxRight, dyRight]
        ]
}

function fn(x) {
    return Math.tan(x);
}

EDIT: Forgot to mention that it's assumed the approximate root was already found within the interval xmin and xmax, and it does exist somewhere between that interval.

So the maximum interval will be the distance from xmin to xmax.

I guess this question is subtle as there is a trade-off between closeness (precision) and iterations (performance), and overall I'm looking for something more optimal than what I have

C9C
  • 319
  • 3
  • 14
  • This is actually a pretty subtle question. Without going into much more detail about floating point numbers, I think a simple and effective approach is to use whatever method you want to get an approximate root, then make a small interval around that and apply bisection. Note that it's not guaranteed that the initial approximate root will be within the final small interval! Good luck and have fun. – Robert Dodier Jul 31 '20 at 21:57
  • @RobertDodier Thanks for the suggestion. I guess that approach makes the most sense for this kind of problem! – C9C Jul 31 '20 at 23:05

1 Answers1

2

I presume that you refer to the floating-point representation of the numbers (because for a continuous function, the "nearest" points do not exist).

Beware that very close to a root the function evaluation loses any reliability and might not even be monotonous and there could be several nearby changes of signs, making the "root" not unique. None of the usual methods, except dichotomy, are even guaranteed to land close to a change of sign !

Anyway, if you want the nearest point of a different sign (chances of an exact zero are pretty low), you can explore on either sides by incrementing/decrementing the lower-order bits of the floating-point mantissa and evaluate the function. The best is to map an integer on the FP representation.

Values that cause a change in the exponent will need special handling (though for first experiments you can probably ignore them).