mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.49k stars 897 forks source link

Brent root finder fails to find roots #1001

Closed naasking closed 1 year ago

naasking commented 1 year ago

Hello, ran into a situation where Brent's root finder failed to find roots, which you can see in the following program:

    static void Main(string[] args)
    {
        var x = new[] {      0,  845.4,   1703,   2459,   3389, 4285.3, 5006.2, 5982.4, 6710.9, 7593.5, 8188.1, 9181.2, };
        var y = new[] { 1.8853, 1.8618, 1.8414, 1.8213, 1.8013, 1.9701, 1.8326, 1.6421,  1.415, 1.0214, 0.6688,      0, };
        var f = CubicSpline.InterpolateNaturalSorted(x, y);
        var y0 = 1.79;
        var x0 = Brent.FindRoot(x1 => f.Interpolate(x1) - y0, 3000, 5800, 0.001, 1000);
        Console.WriteLine($"x0 = {x0}, y0 = {y0}, y(interpolated) = {f.Interpolate(x0)}");
        Console.ReadLine();
    }

There's a clear peak at x=4285.3, y=1.97, but no matter how many iterations I specify or how accurate/inaccurate the match, it can't find a root. It can find a root if I narrow the search range to 3800-5800, or if the y0 < 1.79.

Not sure if this is expected with this method, but it was surprising.

jkalias commented 1 year ago

Hi,

It seems you hit a nasty corner with this choice of points and target value y0. According to the implemenation you need to have a changing sign so that the algorithm works.

If you plot the spline interpolation, you will see that for 1.79 you have multiple roots, and your initial range selection (3000, 5800) hits two y points, which have the same sign. image

I'm not sure this helps you

naasking commented 1 year ago

Interesting, thanks. I didn't see that described in the documentation and hadn't looked into the code. I had already resolved this by switching to Newton-Raphson awhile back, but thought it might be worth lodging a ticket. Feel free to close if this isn't going to influence the code at all.

jkalias commented 1 year ago

I agree that this is a bit frustrating, but it is a common limitation of bracketing methods.

I don’t see a valid way of enhancing a general purposing library like Numerics to deal with such issues. You can safely close the ticket, thanks for the effort and your time.