我正在尝试加快我的代码,目前在 Python/Numpy 中运行它需要一个多小时。大部分计算时间发生在下面粘贴的函数中。
我正在尝试矢量化 Z,但我发现它对于三重 for 循环来说相当困难。我可以在某处实现 numpy.diff
函数吗?看一看:
def MyFESolver(KK,D,r,Z):
global tdim
global xdim
global q1
global q2
for k in range(1,tdim):
for i in range(1,xdim-1):
for j in range (1,xdim-1):
Z[k,i,j]=Z[k-1,i,j]+r*q1*Z[k-1,i,j]*(KK-Z[k-1,i,j])+D*q2*(Z[k-1,i-1,j]-4*Z[k-1,i,j]+Z[k-1,i+1,j]+Z[k-1,i,j-1]+Z[k-1,i,j+1])
return Z
tdim = 75
xdim = 25
最佳答案
我同意,这很棘手,因为所有四个边上的 BC 破坏了刚度矩阵的简单结构。您可以这样摆脱空间循环:
from pylab import *
from scipy.sparse.lil import lil_matrix
tdim = 3; xdim = 4; r = 1.0; q1, q2 = .05, .05; KK= 1.0; D = .5 #random values
Z = ones((tdim, xdim, xdim))
#Iterate in time
for k in range(1,tdim):
Z_prev = Z[k-1,:,:] #may need to flatten
Z_up = Z_prev[1:-1,2:]
Z_down = Z_prev[1:-1,:-2]
Z_left = Z_prev[:-2,1:-1]
Z_right = Z_prev[2:,1:-1]
centre_term = (q1*r*(Z_prev[1:-1,1:-1] + KK) - 4*D*q2)* Z_prev[1:-1,1:-1]
Z[k,1:-1,1:-1]= Z_prev[1:-1,1:-1]+ centre_term + q2*(Z_up+Z_left+Z_right+Z_down)
但我不认为你可以摆脱时间循环...
我认为表达:
Z_up = Z_prev[1:-1,2:]
在 numpy 中制作一个副本,而你想要的是一个 View - 如果你能弄清楚如何做到这一点 - 它应该更快(多少?)
最后,我同意其他回答者的观点——根据经验,这种循环最好用 C 语言完成,然后再封装到 numpy 中。但是上面应该比原来的要快...
关于python - 是否可以在 Python/Numpy 中对这个三重 for 循环进行矢量化?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13065596/