2013-04-14 2 views
5
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 

이 코드의 잘못된 점은 무엇입니까? 그것은 작동하지 않는 것 같습니다.솔루션을 찾으려면 fsolve를 사용하십시오.

+2

"작동하지 않음"을 정의하십시오. –

+0

두 번째'sqrt' 매개 변수의 분모를 지정하는 데 오류가있는 것 같습니다. 아마도'np.sqrt (.5 * musun ** 2/(e * (1 - x ** 2))))'? – mtadd

답변

4

fsolve()f(x) = 0 (here 참조)의 뿌리를 반환합니다. I의 범위에 대한 xf(x)의 값을 플롯

-1 (1), I는 x = -1x = 1에서 뿌리가 있다는 것을 발견했다. 그러나 x > 1 또는 x < -1 인 경우 sqrt() 함수 모두에 음수 인수가 전달되어 오류 invalid value encountered in sqrt이 발생합니다.

fsolve()은 함수의 유효 범위의 맨 끝에있는 루트를 찾는 데 실패하지 않습니다.

나는 그 뿌리를 찾으려고하기 전에 함수의 그래프를 그려 보는 것이 좋다는 것을 알 수있다. 그 이유는 그것이 그 뿌리를 찾을 확률이 얼마나 높을 지 (또는이 경우에는있을 법하지 않음) 나타낼 수 있기 때문이다. 모든 루트 찾기 알고리즘.

8

sqrt은 nagative 인수로 NaN을 반환하므로 함수 f (x)는 모든 실제 x에 대해 계산할 수 없습니다. 나는 < 인수가 0 일 때 복소수 값을 출력 할 수있는 함수 numpy.emath.sqrt()을 사용하여 표현식의 절대 값을 반환합니다.

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.

관련 문제