用于二维热方程的 numpy 快速方法答案

【问题标题】：fast method with numpy for 2D Heat equation用于二维热方程的 numpy 快速方法
【发布时间】：2015-05-30 18:15:39
【问题描述】：

我正在寻找一种用 python 求解二维热方程的方法。我已经实现了有限差分法，但是是慢动作（进行 100,000 次模拟需要 30 分钟）。这个想法是创建一个端可以写的代码，

for t in TIME:  
    DeltaU=f(U)
    U=U+DeltaU*DeltaT
    save(U)

我该怎么做？

在我的代码的第一种形式中，我使用了有限差分的 2D 方法，我的格栅是 5000x250 (x, y)。现在我想降低计算速度，想法是找到

DeltaU = f(u)

其中 U 是热函数。为了实现，我将此源 http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/ 用于 2D 案例，但运行时间因我的需要而更加昂贵。有没有一些方法可以做到这一点？

也许我必须使用矩阵

A=1/dx^2 (2 -1  0  0 ... 0 
         -1  2 -1  0 ... 0
         0  -1  2 -1 ... 0
         .            .
         .            .
         .            .
         0 ...        -1 2)

但是如何在 2D 问题中做到这一点？如何在A中插入边界条件？这是我使用的有限差分代码：

Lx=5000   # physical length x vector in micron
Ly=250   # physical length y vector in micron
Nx = 100     # number of point of mesh along x direction
Ny = 50     # number of point of mesh along y direction
a = 0.001 # diffusion coefficent
dx = 1/Nx
dy = 1/Ny
dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2)) # it is 0.04
x = linspace(0.1,Lx, Nx)[np.newaxis] # vector to create mesh
y = linspace(0.1,Ly, Ny)[np.newaxis] # vector to create mesh
I=sqrt(x*y.T) #initial data for heat equation
u=np.ones(([Nx,Ny])) # u is the matrix referred to heat function
steps=100000
for m in range (0,steps):
        du=np.zeros(([Nx,Ny]))

        for i in range (1,Nx-1):

            for j in range(1,Ny-1):

               dux = ( u[i+1,j] - 2*u[i,j] + u[i-1, j] ) / dx**2
               duy = ( u[i,j+1] - 2*u[i,j] + u[i, j-1] ) / dy**2            
               du[i,j] = dt*a*(dux+duy)




    # Boundary Conditions
    t1=(u[:,0]+u[:,1])/2
    u[:,0]=t1
    u[:,1]=t1
    t2=(u[0,:]+u[1,:])/2
    u[0,:]=t2
    u[1,:]=t2
    t3=(u[-1,:]+u[-2,:])/2
    u[-1,:]=t3
    u[-2,:]=t3
    u[:,-1]=1


    filename1='data_{:08d}.txt'

    if m%100==0:
        np.savetxt(filename1.format(m),u,delimiter='\t' )

对于复杂的 100000 步，运行时间约为 30 分钟。我将优化此代码（使用初始行中提出的想法）以使运行时间约为 5/10 分钟或更短。我该怎么做？

【问题讨论】：

欢迎来到堆栈溢出！您的问题就目前而言非常广泛，而且很难理解具体问题在哪里。你也问过an identical question，两次问同一个问题通常是个坏主意。您在 cmets 中就该问题提供了一些很好的建议，除非您更详细地展示您的确切实现以及您认为它作为minimal, complete and verifiable example 太慢的地方，否则您不太可能获得更好的帮助
对不起，我是新手，不知道这个论坛的指导方针。我只是想知道在这个想法的基础上是否有比有限差分更快的方法以上概述。如果我提出了两个相同的问题，我深表歉意，但我无法理解第一个申请的暂停。现在我该怎么处理这个问题？
没问题，感谢您的回复。我认为您可以将其编辑为更好的问题以寻求帮助。包括实际执行有限差分的代码部分、您计算的点数（即您的网格大小）以及它的运行速度与您认为它可以/希望它运行的速度有多快
然后，打开另一个问题或对此发表评论？
编辑这个（你仍然可以在它被搁置时编辑）。这是一个建议。显示执行循环和有限差分逻辑的确切代码。您的版本，而不是简单的链接，将代码放入问题中。查找 python timeit 模块并根据您的情况发布结果。写下为什么您认为基于矩阵的解决方案可能会更快以及您遇到的问题（即为什么您不能直接去尝试）。那么它可能有资格重新开放

标签： python python-3.x numpy heat

【解决方案1】：

您是否考虑过并行化您的代码或使用 GPU 加速。

如果您在 python 分析器 (cProfile) 中运行代码会有所帮助，这样您就可以找出运行时的瓶颈所在。我假设它是在求解你得到的矩阵方程，可以通过我上面列出的方法轻松加速。

【讨论】：

代码工作正常，我必须让它更快，我认为使用函数 DeltaU=f(u) 来减少运行时间。此外，代码应该在慢速机器上运行，因此我想对其进行优化。

【解决方案2】：

有一些简单但巨大的改进可能：

仅通过引入“Dxx, Dyy = 1/(dxdx), 1/(dydy)”，运行时间下降 25%。通过使用切片并避免 for 循环，代码现在 快了 400 倍

def f_for(u):
    for m in range (0,steps):
            du=np.zeros((Nx,Ny))
            for i in range (1,Nx-1):
                for j in range(1,Ny-1):
                   dux = ( u[i+1,j] - 2*u[i,j] + u[i-1, j] ) / dx**2
                   duy = ( u[i,j+1] - 2*u[i,j] + u[i, j-1] ) / dy**2            
                   du[i,j] = dt*a*(dux+duy)
                
                
def f_slice(u):
    Dxx, Dyy = 1/dx**2, 1/dy**2
    i = slice(1,nx-1); iw = slice(0,nx-2); ie = slice(2,nx)
    j = slice(1,ny-1); js = slice(0,ny-2); jn = slice(2,ny)
    for m in range (0,steps):
        dux = Dxx* ( u[ie,j] - 2*u[i,j] + u[iw, j] )
        duy = Dyy* ( u[i,jn] - 2*u[i,j] + u[i, js] )           
        du[i,j] = dt*a*(dux+duy)
            
            
Nx = 100     # number of point of mesh along x direction
Ny = 50     # number of point of mesh along y direction
a = 0.001 # diffusion coefficent
dx = 1/Nx
dy = 1/Ny
dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2))
steps=100000
steps=100
U=np.ones(([Nx,Ny])) # u is the matrix referred to heat function

%timeit f_for(U)
%timeit f_slice(U)

【讨论】：

【解决方案3】：

我可能错了，但在您为时间步创建的循环中的代码中，“m in range(steps)” 在下面的一行中继续；杜=np.zeros(----)。不是 python 专家，但这可能会导致创建一个稀疏矩阵的步数，在本例中为 100k 次。

【讨论】：