使用fsolve找到解决方案

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

这段代码有什么问题？它似乎行不通.

Answer 1

因为sqrt为nagative参数返回NaN,所以函数f(x)不能计算所有实数x.我将你的函数改为使用numpy.emath.sqrt(),当参数<0时,它可以输出复数值,并返回表达式的绝对值.

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)

然后你就可以得到正确的结果:

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

解决方案非常接近真正的根,但f(x)仍然非常大,因为f(x)有一个非常大的因素:musun.

Answer 2

fsolve()返回f(x) = 0(见这里)的根.

当我绘制的值f(x)对x在范围-1到1,我发现有在根x = -1和x = 1.但是,如果x > 1或者x < -1,两个sqrt()函数都将传递一个负参数,这会导致错误invalid value encountered in sqrt.

我没有惊讶的是fsolve(),找不到函数有效范围最末端的根.

我发现在尝试查找函数的根之前绘制函数的图形总是一个好主意,因为这可以表明任何根查找发现根的可能性(或者在这种情况下,不太可能)算法.