1

題名の通り、\vec{z}'=fという2階の微分方程式(数値的に解くために一階のベクトルの微分方程式の形にしている)を、初期値X0とパラメータepsilonを変化させて、無限遠で0になるという境界条件に合う(x0,epsilon)の組み合わせを見つけて、その時のepsilonの最小値を拾ってくるということをpythonでしたいのですが、以下のようなコードを書いたところ、やはりループがいけないのかとっても時間がかかってしまいます。本当はmeshgridなどで出来れば一番早いと思うのですが、うまくいきませんでした。

極力ループを使わないような方法で以下のコードと同じようなことができる方法があればご教授お願い致します。

def sol(q_x,kappa,h,x0,epsilon):

    y = Symbol('y')

    def fsol(x,y,epsilon):
        return f(x,y,q_x,kappa,h,epsilon)

    y = np.linspace(0.0,1000,10001)

    return odeint(fsol,x0,y,args=(epsilon,))

alpha = np.linspace(0,2*np.pi,101)
X = np.cos(alpha)
Y = np.sin(alpha)
ZEROS = [[0.0]*2 for i in range(len(alpha))]
X0 = np.c_[X,Y,ZEROS]

def dispersion0(q_x,kappa,h,cutoff):

    def z(x0,epsilon):
        return sol(q_x,kappa,h,x0,epsilon)

    epsilon = np.linspace(0,q_x**2+kappa+h,101)

    Phi1 = np.empty([alpha.shape[0],epsilon.shape[0]])
    Phi2 = np.empty([alpha.shape[0],epsilon.shape[0]])
    for i in range(0,len(alpha)):
        for j in range(0,len(epsilon)):
            Phi1[i][j] = z(X0[i],epsilon[j])[1000][0]
            Phi2[i][j] = z(X0[i],epsilon[j])[1000][1]
    energy = 0
    for k in range(0,len(epsilon)):
        if abs(Phi1[:,k])<cutoff & abs(Phi2[:,k])<cutoff :
            energy = epsilon(k)
            break            

    return energy

1 件の回答 1

1

質問のコードをみると、solは、常微分方程式を解くルーティンですが、その関数を、dispersion0から、2百万回呼んでいます。solの計算に1msかかったとしても、33分必要なのですが、vec{z}についての具体的な説明がないのでハッキリしないところもありますが、PCでsolの計算が1msでは終わるとは通常考えられないので、このままでは処理は終わらないと思います。

解決策としてはdispersion0で総当りで計算していますが、それを減らすしかないと思われます。数値解析の反復解法を使う、探索木での枝刈りを使って解になりそうにないものを除いてしまうという方法があります。

また、fsolは呼び出される回数が非常に多くなるので、SimPyのような処理の遅いものを使わずにCythonを使って高速化すべきです。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。