using fsolve to find the solution

Question

import numpy as np
from scipy.optimize import fsolve
musun = 132712000000
T = 365.25 * 86400 * 2 / 3
e = 581.2392124070273
def f(x):
 return ((T * musun ** 2 / (2 * np.pi)) ** (1 / 3) * np.sqrt(1 - x ** 2)
 - np.sqrt(.5 * musun ** 2 / e * (1 - x ** 2)))
x = fsolve(f, 0.01)
f(x)
print x

What is wrong with this code? It seems to not work.

• mtadd: It looks like there might be an error in specifying the denominator of your second sqrt parameter. Perhaps np.sqrt(.5 * musun ** 2 / (e * (1 - x ** 2))))?

Answer 1

Because sqrt returns NaN for negative argument, your function f(x) is not calculable for all real x. I changed your function to use numpy.emath.sqrt() which can output complex values when the argument < 0, and returns the absolute value of the expression.

import numpy as np
from scipy.optimize import fsolve
sqrt = np.emath.sqrt
musun = 132712000000
T = 365.25 * 86400 * 2 / 3
e = 581.2392124070273
def f(x):
 return np.abs((T * musun ** 2 / (2 * np.pi)) ** (1 / 3) * sqrt(1 - x ** 2)
 - sqrt(.5 * musun ** 2 / e * (1 - x ** 2)))
x = fsolve(f, 0.01)
x, f(x)

Then you can get the right result:

(array([ 1.]), array([ 121341.22302275]))

the solution is very close to the true root, but f(x) is still very large because f(x) has a very large factor: musun.

Answer 2

fsolve() returns the roots of f(x) = 0 (see here).

When I plotted the values of f(x) for x in the range -1 to 1, I found that there are roots at x = -1 and x = 1. However, if x > 1 or x < -1, both of the sqrt() functions will be passed a negative argument, which causes the error invalid value encountered in sqrt.

It doesn't surprise me that fsolve() fails to find roots that are at the very ends of the valid range for the function.

I find that it is always a good idea to plot the graph of a function before trying to find its roots, as that can indicate how likely (or in this case, unlikely) it is that the roots will be found by any root-finding algorithm.