題名の通り、\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