CODEWITHSUNDEEP
pythonoptimizationscipypolynomial-math
Alex
Alex

Reputation: 4264

Solve for roots in given interval using scipy.optimize

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?

Upvotes: 0

Views: 5053

Answers (3)

Bill Bell
Bill Bell

Reputation: 21663

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.

Function to study

Here's some typical code in this case.

Produce plot and find roots

Upvotes: 1

Warren Weckesser
Warren Weckesser

Reputation: 114921

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

Upvotes: 2

ev-br
ev-br

Reputation: 26080

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.

Upvotes: 2

Related Questions