

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

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

[code]

    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
    
[/code]


CC BY-SA 4.0
2018-09-10T08:20:01.977
48255
2020-08-26T09:06:10.780
2019-06-21T02:05:48.020
3060
30055
question
python
初期値と固有値を変化させて微分方程式を解く(Shooting method)のpythonでの実装について


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

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

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


CC BY-SA 4.0
2018-09-12T05:29:19.393
48304
2018-09-12T05:29:19.393
15171
48255
answer