2017-02-22 4 views
2

n(ポリトロープインデックス)の任意の値に対するLane-Emden方程式を解こうとしています。 SciPyを使用するために、2次ODEを2つの結合された1次ODEの集合として表現しました。私は、まずSciPyでODEを解くとdouble_scalarsエラーで無効な値が返される

RuntimeWarning: invalid value encountered in double_scalars
app.launch_new_instance()

:私は次のエラーで、このコードの結果を実行している、しかし、ファイ=シータ」

y0 = [1.,0.] 
xi = np.linspace(0., 16., 201)  
for n in range(0,11): 
    sol = odeint(poly, y0, xi, args=(n/2.,)) 
    plt.plot(xi, sol[:, 0], label=str(n/2.)) 
plt.legend(loc='best') 
plt.xlabel('xi') 
plt.grid() 
plt.show() 

を定義し、ここで

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 

def poly(y,xi,n): 
    theta, phi = y 
    dydt = [phi/(xi**2), -(xi**2)*theta**n] 
    return dydt 

:私は、次のコードを持っていますこれはxi = 0での方程式の特異点の結果だと思ったので、積分間隔を変更しました:

xi = np.linspace(1e-10, 16., 201) 

これはn = 0の問題を修正しますが、nの他の値は修正しません。私が得たプロットは意味をなさないだけで間違っています。

なぜこのエラーメッセージが表示され、コードを修正できますか?

答えて

1

以下は私のために働きます(Wikipedia entryと似ています)。その派生語は誤って書かれていた(間違った要素に負である)、派生語への入力は単にnに変更された。

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import odeint 

def poly(y,xi,n): 
    theta, phi = y 
    dydxi = [-phi/(xi**2), (xi**2)*theta**n] 
    return dydxi 

fig, ax = plt.subplots() 

y0 = [1.,0.] 
xi = np.linspace(10e-4, 16., 201) 

for n in range(6): 
    sol = odeint(poly, y0, xi, args=(n,)) 
    ax.plot(xi, sol[:, 0]) 

ax.axhline(y = 0.0, linestyle = '--', color = 'k') 
ax.set_xlim([0, 10]) 
ax.set_ylim([-0.5, 1.0]) 
ax.set_xlabel(r"$\xi$") 
ax.set_ylabel(r"$\theta$") 
plt.show() 

​​


元の質問ポリトロープ指数の端数値を使用するので、フォローのようなものは、それらの値を調査するために使用することができ...

for n in range(12): 
    sol = odeint(poly, y0, xi, args=(n/2.,)) 

    if np.isnan(sol[:,1]).any(): 
     stopper = np.where(np.isnan(sol[:,1]))[0][0] 
     ax.plot(xi[:stopper], sol[:stopper, 0]) 
    else: 
     ax.plot(xi, sol[:, 0]) 

fractional polytropic index values

+0

返信ありがとうございます!それは素晴らしいです、私のコードは今動作します!なぜそれがnの整数値に対してしか機能しないのか分かりますか? –

+1

'theta 'がゼロ以下のときは好きではないと思います。解の出力を見ると、それは 'theta = 0'で停止します。 –

+1

これは理にかなっています。複雑な解決策につながるからです。 theta <0のポイントまでODEintを使用する方法はありますか?興味深いもの: –

関連する問題