6
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

这段代码有什么问题?它似乎不起作用。

4

2 回答 2

12

因为sqrt为负参数返回 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。

于 2013-04-14T12:24:48.867 回答
7

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

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

fsolve()未能找到位于函数有效范围末端的根并不令我感到惊讶。

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

于 2013-04-14T06:06:49.190 回答