używanie fsolve do znalezienia rozwiązania

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 

Co jest nie tak z tym kodem? Wygląda na to, że nie działa.

Answer 1

fsolve() zwraca pierwiastki z f(x) = 0 (patrz here).

Kiedy wykreślono wartości f(x) do x w zakresie od -1 do 1, okazało się, że nie są korzenie w x = -1 i x = 1. Jednakże, jeśli x > 1 lub x < -1, zarówno funkcji sqrt() zostanie przekazany argument ujemny, co powoduje błąd invalid value encountered in sqrt.

To mnie nie dziwi, że fsolve() nie znajdzie korzenie, które znajdują się na krańce prawidłowego zakresu dla funkcji.

Uważam, że zawsze dobrym pomysłem jest wykreślenie wykresu funkcji przed próbą znalezienia jej korzeni, ponieważ może to wskazywać na prawdopodobieństwo (lub, w tym przypadku, nieprawdopodobne), że korzenie zostaną znalezione przez dowolny algorytm odnajdywania rootów.

Answer 2

Ponieważ sqrt zwraca NaN dla nagative argumentu, to funkcja f (x) nie jest calculatable dla wszystkich rzeczywistych x. Zmienić swoją funkcję do korzystania numpy.emath.sqrt() które mogą wyjściowych wartości zespolone, gdy argument < 0 i zwraca wartość bezwzględną wyrażenia.

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) 

Następnie można uzyskać odpowiedni wynik:

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

rozwiązanie jest bardzo blisko do prawdziwego korzenia, ale f (x) jest nadal bardzo duża, ponieważ f (x) ma bardzo duży czynnik: musun.