используя 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

fsolve() возвращает корни f(x) = 0 (см. here).

Когда я отложены значения f(x) для x в диапазоне от -1 до 1, я обнаружил, что есть корни в x = -1 и x = 1. Однако, если x > 1 или x < -1, обе функции sqrt() будут переданы отрицательным аргументом, что вызывает ошибку invalid value encountered in sqrt.

Меня не удивляет, что fsolve() не находит корни, которые находятся в самом конце допустимого диапазона для функции.

Я считаю, что всегда полезно построить график функции, прежде чем пытаться найти ее корни, поскольку это может указывать, насколько вероятно (или в этом случае маловероятно), что корни будут найдены любой алгоритм поиска корней.

Answer 2

Поскольку sqrt возвращает NaN для аргумента nagative, функция 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 (х) имеет очень большой фактор: musun.