解決法を見つけるために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を参照）。 Iが範囲内xためf(x)の値をプロット

は

-1 1に、Iはx = -1とx = 1で根が存在することを見出しました。ただし、x > 1またはx < -1の場合、sqrt()関数の両方に負の引数が渡されます。これにより、エラーinvalid value encountered in sqrtが発生します。

fsolve()は、関数の有効範囲の最後にあるルートを見つけることができません。

私は、その根を見つける前に、関数のグラフをプロットすることをお勧めします。それは、根が見つかる可能性が高い（この場合はそうではない）ことを示すことができるからです。任意のルート探索アルゴリズム。

Answer 2

sqrtはnagative引数として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。