utilizzando fsolve per trovare la soluzione

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 

Cosa c'è di sbagliato con questo codice? Sembra non funzionare.

Answer 1

fsolve() restituisce le origini di f(x) = 0 (vedere here).

Quando ho tracciato i valori di f(x) per x nell'intervallo da -1 a 1, ho trovato che ci sono radici alla x = -1 e x = 1. Tuttavia, se x > 1 o x < -1, a entrambe le funzioni sqrt() verrà passato un argomento negativo, che provoca l'errore invalid value encountered in sqrt.

Non mi sorprende che fsolve() non riesca a trovare le radici che si trovano all'estremità dell'intervallo valido per la funzione.

Trovo che sia sempre una buona idea tracciare il grafico di una funzione prima di provare a trovare le sue radici, in quanto ciò può indicare quanto sia probabile (o in questo caso improbabile) che le radici vengano trovate da qualsiasi algoritmo di ricerca delle radici.

Answer 2

Perché sqrt restituisce NaN per argomento nagativo, la funzione f (x) non è calcolabile per tutto il x reale. Cambio la tua funzione per usare numpy.emath.sqrt() che può emettere valori complessi quando l'argomento < 0 restituisce il valore assoluto dell'espressione.

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) 

allora si può ottenere il risultato giusto:

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

la soluzione è molto vicino alla vera radice, ma f (x) è ancora molto grande, perché f (x) ha un grande fattore: musun.