Solve for roots in given interval using scipy.optimize

Question

I have the function f1 = lambda x: 1 - 1.12 * (x ** 0.5) * ((1-x) ** 0.02), and I wish to solve for its roots in the interval (0, 1). I have tried using the scipy.optimize.brentq and scipy.optimize.fsolve to do this, but both methods run into issues. Based on some experimentation, I got that the roots of this equation are approximately equal to 0.86322414 and 0.9961936895432034 (we know there are at most 2 roots because the function has one inflection point in this interval):

f1 = lambda x: 1 - 1.12 * (x ** 0.5) * ((1-x) ** 0.02)
print(fsolve(f1, 0.5))
print(f1(0.99))
print(f1(0.999))
print(brentq(f1, 0.99, 0.999))

Output:

[ 0.86322414]
-0.016332046983897452
0.025008640855473052
0.9961936895432034

The issue here is that in order for brentq to work, the values of the function must be of opposite signs at the specified endpoints. Furthermore, when I started fsolve at values of x close to 1, I got runtime warning messages:

print(fsolve(f1, 0.97))
print(fsolve(f1, 0.98))

Output:

[ 0.97]
[ 0.98]
C:/Users/Alexander/Google Drive/Programming/Projects/Root Finding/roots.py:6: RuntimeWarning: invalid value encountered in power

C:\Users\Alexander\Anaconda3\lib\site-packages\scipy\optimize\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.

Does anyone if there is a more systematic way to solve for roots of this equation, and why fsolve is not working on x = 0.97, 0.98?

FWIW, using mpmath's `findroot` with the 'mnewton' solver I get 0.863224139030458644097939619031454833520588010 and 0.996193689543209644983205954153079782870170668 — PM 2Ring, Aug 27 '16 at 15:40
Here is one solution which might help you http://stackoverflow.com/questions/28187569/find-the-root-of-a-cubic-function — Abhijay Ghildyal, Aug 27 '16 at 15:42
@PM2Ring can you go into a bit more detail about how you used findroot with the mnewton solver? I tried findroot(f1, 0.5, solver = "mnewton") but that only gave me one solution. Thanks. — Alex, Aug 27 '16 at 20:49
I got the higher root with `mp.findroot(f1, 0.97, solver='mnewton')`. I admit that you have to use trial & error to find a starting approximation that will lead to convergence on that root, and using 1 as the initial approximation raises `ValueError: Could not find root within given tolerance.` Finding roots where there isn't a sign change is tricky! BTW, when using mpmath at high precision it's important to create `mpf` numbers via strings or integers rather than Python floats. See [Providing correct input](http://mpmath.org/doc/current/basics.html#providing-correct-input) — PM 2Ring, Aug 28 '16 at 08:02

score 2 · Answer 1 · edited Aug 28 '16 at 02:16

fsolve does not support root-finding on an interval. It likely wanders off to either x<0 or x>1, hence RuntimeWarning (once) and garbage for the answer. You can check it by instrumenting your function with print(x).

Brentq's behavior is undefined if the interval contains more then one root.

If you know the inflection point, or know that there's only one, then you can find it via brentq and use it to bracket your two roots.

score 2 · Accepted Answer · answered Aug 28 '16 at 04:47

If you take the derivative of your function and set it equal to 0, after a little algebra you'll find that the derivative is 0 at x0 = 0.5/0.52. (In a calculus class, this point is called a critical point, not an inflection point.) The function has a minimum at this point, and the value there is negative. The values at x=0 and x=1 are positive, so you can use [0, x0] and [x0, 1] as bracketing intervals in brentq:

In [17]: from scipy.optimize import brentq

In [18]: f1 = lambda x: 1 - 1.12 * (x ** 0.5) * ((1-x) ** 0.02)

In [19]: x0 = 0.5/0.52

In [20]: brentq(f1, 0, x0)
Out[20]: 0.8632241390303161

In [21]: brentq(f1, x0, 1)
Out[21]: 0.9961936895432096

score 1 · Answer 3 · answered Aug 28 '16 at 19:01

As you already know, a direct answer to part of your question is that fsolve doesn't find a root in the interval [0.97,0.98] because there's no root there. As for a systematic way, why not use plot?

Once you've defined the function as a lambda, ready for use with Brent or various other routines, you can just copy and paste the substring you need into a call to plot and indicate the interval of interest. If one of the roots is a bit obscure then broaden that part of the x-axis.

Here's some typical code in this case.

Solve for roots in given interval using scipy.optimize

3 Answers3