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.

Answer 1

Because sqrt returns NaN for nagative argument, you function f(x) is not calculatable for all real x. I change 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.