python - 用python对电势进行数值微分

Question

我使用以下代码获得四极电位的电场。

for n in range (nx-1):
  for m in range (ny-1):
   for k in range (nz-1):
     Ex[n,m,k] =-( (T[n+1,m,k]-T[n,m,k]) )/((x[n+1]-x[n]));
     Ey[n,m,k] =-( (T[n,m+1,k]-T[n,m,k]) )/((y[m+1]-y[m]));
     Ez[n,m,k] =-( (T[n,m,k+1]-T[n,m,k]) )/((z[k+1]-z[k]));
     return Ex,Ey,Ez, T

这里 T 是通过数值求解拉普拉斯方程获得的 3D 电势，从图中可以看出，红色电极(正极)的电场矢量方向错误(顶部电极的右上侧和右下侧)，其他电极也有同样的错误。即负电极有出射电场矢量，方向必须相反。

我也使用了中心差分法，但我得到了相同的数字。你能告诉我我的差异化有什么问题吗？

Answer 1

您的问题在于 matplotlib 解释方向的方式，这本身源于矩阵的常规索引方式，这与许多人(例如我自己)的想法相反。具体来说，第一个索引是垂直 索引，而第二个索引是水平 索引。如果您使用的是非方形数组，这会很清楚，这会显示您的 x 和 y 是向后的。以下示例说明了您的渐变(仅稍作修改)如何给出正确的结果。第一个图显示了您实际绘制的内容，其中交换了渐变的 x 和 y 分量。该图(您可以通过运行它获得)表明渐变与等高线不正交(因为它必须如此)并且有时会朝错误的方向移动。

第二个和第三个图都显示了使用很好的方法绘制梯度的正确方法，在该方法中我们有坐标的二维数组以及我们正在绘制的东西(无论是 quiver 还是轮廓)。我按照您的代码使用 meshgrid 生成我的 X 和 Y，但需要将参数交换到 meshgrid获得正确的尺寸。一种更简洁的方法是使用相同的编码风格来生成与绘制的事物相同的坐标，在这种情况下，具有显式索引的 while 循环将是最佳的。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-3,3,0.2)
y = np.arange(-5,5,0.2)
z = np.arange(-5,5,0.2)

def grad(f):
    gx = np.zeros_like(f)
    gy = np.zeros_like(f)
    gz = np.zeros_like(f)
    for n in range(1,len(x)-1):
        for m in range(1,len(y)-1):
            for k in range(1,len(z)-1):
                gx[n,m,k] = (T[n+1,m,k]-T[n-1,m,k])/(x[n+1]-x[n-1]);
                gy[n,m,k] = (T[n,m+1,k]-T[n,m-1,k])/(y[m+1]-y[m-1]);
                gz[n,m,k] = (T[n,m,k+1]-T[n,m,k-1])/(z[k+1]-z[k-1]);
    return gx, gy, gz

T = np.zeros((len(x), len(y), len(z)))
for n in range(len(x)):
    for m in range(len(y)):
        for k in range(len(z)):
            T[n,m,k] = np.sin((x[n] - y[m])/3.0) + 0.3*np.cos(y[m]) + z[k]**2

gx,gy,gz = grad(T)

Y, X= np.meshgrid(y,x)

plt.figure('WRONG with x and y')
plt.contour(y, x, T[:,:,round(len(z)/2)], 64)
plt.colorbar()
plt.quiver(y, x, 10*gx[:,:,round(len(z)/2)], 10*gy[:,:,round(len(z)/2)])
plt.xlabel("x")
plt.ylabel("y")
plt.axes().set_aspect('equal')


plt.figure('with X and Y')
plt.contour(X, Y, T[:,:,round(len(z)/2)], 64)
plt.colorbar()
plt.quiver(X, Y, 10*gx[:,:,round(len(z)/2)], 10*gy[:,:,round(len(z)/2)])
plt.xlabel("X")
plt.ylabel("Y")
plt.axes().set_aspect('equal')

plt.figure('with Y and X')
plt.contour(Y, X, T[:,:,round(len(z)/2)], 64)
plt.colorbar()
plt.quiver(Y, X, 10*gy[:,:,round(len(z)/2)], 10*gx[:,:,round(len(z)/2)])
plt.xlabel("Y")
plt.ylabel("X")
plt.axes().set_aspect('equal')

plt.show()

最后，我要再次指出，揭示此错误的最简洁方法是使用 nx != ny 运行您的程序。该代码会失败，并显示一条错误消息，指示数组维度不匹配，这会导致您(希望)以正确的方式交换事物。