PDE の有限差分法を表す行列がAx=bどこにあるかを解決する必要があります。A2D 問題の典型的なサイズはA(256^2)x(256^2) 前後で、いくつかの対角線で構成されています。次のサンプル コードは、A の作成方法です。
N = Nx*Ny # Nx is no. of cols (size in x-direction), Ny is no. rows (size in y-direction)
# finite difference in x-direction
up1 = (0.5)*c
up1[Nx-1::Nx] = 0
down1 = (-0.5)*c
down1[::Nx] = 0
matX = diags([down1[1:], up1[:-1]], [-1,1], format='csc')
# finite difference in y-direction
up1 = (0.5)*c
down1 = (-0.5)*c
matY = diags([down1[Nx:], up1[:N-Nx]], [-Nx,Nx], format='csc')
とを足すmatXとmatY4 つの対角線になります。上記は二次離散化の場合です。4 次離散化の場合、8 つの対角線があります。二次導関数がある場合、主対角もゼロではありません。
次のコードを使用して、実際の解決を行います。
# Initialize A_fixed, B_fixed
if const is True: # the potential term V(x) is time-independent
A = A_fixed + sp.sparse.diags(V_func(x))
B = B_fixed + sp.sparse.diags(V_func(x))
A_factored = sp.sparse.linalg.factorized(A)
for idx, t in enumerate(t_steps[1:],1):
# Solve Ax=b=Bu
if const in False: #
A = A_fixed + sp.sparse.diags(V_func(x,t))
B = B_fixed + sp.sparse.diags(V_func(x,t))
psi_n = B.dot(psi_old)
if const is True:
psi_new = A_factored(psi_n)
else:
psi_new = sp.sparse.linalg.spsolve(A,psi_n,use_umfpack=False)
psi_old = psi_new.copy()
いくつか質問があります。
Ax=bscipyで解決する最良の方法は何ですか? LU分解を使用するライブラリで使用しますspsolve。sp.sparse.linalg密行列の標準sp.linalg.solveを使用してみましたが、かなり遅くなります。ライブラリ内の他のすべての反復ソルバーも使用してみましたsp.sparse.linalgが、それらも遅くなりました (1000x1000 の場合はすべて数秒かかりますがspsolve、A. 計算を行う別の方法はありますか?
編集:私が解決しようとしている問題は、実際には時間依存のシュレディンガー方程式です。潜在的な項が時間に依存しない場合、提案されているように、最初に行列を事前因数分解Aしてコードを高速化できますが、両方の対角項を変更する必要があるため、潜在的な変数が時変である場合、これは機能しません。行列AとB各時間ステップで。この特定の問題について、事前因数分解または他の方法と同様の方法を使用してコードを高速化する方法はありますか?
- をインストールし
umfpackました。True と False を設定use_umfpackしてテストしてみましたが、驚くべきことに のuse_umfpack=True2 倍近く時間がかかりますuse_umfpack=False。このパッケージを使用すると、速度が向上するはずだと思いました。それが何であるか、何か考えはありますか？(PS: Anaconda Spyder を使用してコードを実行していますが、違いがある場合)
私は自分のコードをプロファイリングするために使用cProfileしており、ほぼすべての時間をspsolve. したがって、コードの残りの部分 (行列/問題の初期化) はかなり最適化されていると思います。改善が必要なのは行列解決の部分です。
ありがとう。