在我们开始之前需要注意一下。range(0, numy-1)
,等于range(numy-1)
,生成从 0 到 numy-2 的数字,不包括 numy-1。那是因为你有从 0 到 numy-2 的 numy-1 值。虽然 MATLAB 具有基于 1 的索引,但 Python 具有基于 0 的索引,因此在转换中的索引时要小心一些。考虑到你有tri_area = np.zeros((numx, numy), dtype=float)
, tri_area[ii,jj]
永远不会以您设置循环的方式访问最后一行或最后一列。因此,我怀疑正确的意图是写range(numy)
.
由于函数t_area()
是可矢量化的,您可以完全消除循环。矢量化意味着 numpy 通过处理底层循环,同时对整个数组应用一些操作,这样它们会更快。
首先,我们将所有的ap
s 表示 (m, n, 3) 数组中的每个 (i, j) 元素,其中 (m, n) 是x
。如果我们计算两个 (m, n, 3) 数组的叉积,则默认情况下该运算将应用于最后一个轴。这意味着np.cross(a, b)
会做对于每个元素 (i, j) 取 3 个数字的叉积a[i,j]
and b[i,j]
。相似地,np.linalg.norm(a, axis=2)
会做对于每个元素 (i, j) 计算 3 个数字的范数a[i,j]
。这也将有效地将我们的数组大小减小到 (m, n)。不过这里要小心一点,因为我们需要明确声明我们希望在第二个轴上完成此操作。
请注意,在以下示例中,我的索引关系可能与您的不对应。完成这项工作的最低要求是surface
多出一行和一列x
and y
.
import numpy as np
def _t_area(a, b, c):
ab = b - a
ac = c - a
return 0.5 * np.linalg.norm(np.cross(ab, ac), axis=2)
def t_area(x, y, surface, dx):
a = np.zeros((x.shape[0], y.shape[0], 3), dtype=float)
b = np.zeros_like(a)
c = np.zeros_like(a)
d = np.zeros_like(a)
a[...,0] = x
a[...,1] = y
a[...,2] = surface[:-1,:-1]
b[...,0] = x + dx
b[...,1] = y
b[...,2] = surface[1:,:-1]
c[...,0] = x
c[...,1] = y + dx
c[...,2] = surface[:-1,1:]
d[...,0] = bp[...,0]
d[...,1] = cp[...,1]
d[...,2] = surface[1:,1:]
# are you sure you didn't mean 0.25???
return 0.5 * (_t_area(a, b, c) + _t_area(d, b, c) + _t_area(b, a, d) + _t_area(c, a, d))
nx, ny = 250, 250
dx = np.random.random()
x = np.random.random((nx, ny))
y = np.random.random((nx, ny))
surface = np.random.random((nx+1, ny+1))
tri_area = t_area(x, y, surface, dx)
x
在此示例中支持索引 0-249,而surface
0-250. surface[:-1]
,简写为surface[0:-1]
,将返回从 0 开始到最后一行的所有行,但不包括它。-1
具有相同的功能并且end
在 MATLAB 中。所以,surface[:-1]
将返回索引 0-249 的行。相似地,surface[1:]
将返回索引 1-250 的行,这与您的效果相同surface[ii+1]
.
Note:我在知道这一点之前就写了这一部分t_area()
可以完全矢量化。因此,虽然这里的内容对于这个答案来说已经过时了,但我将把它作为遗产保留下来,以展示如果函数不可矢量化,可以进行哪些优化。
您应该传递它,而不是为每个元素调用昂贵的函数x
, y,
, surface
and dx
并进行内部迭代。这意味着只需一次函数调用并且开销更少。
此外,您不应该为以下内容创建数组ap
, bp
, cp
and dp
每个循环,这又增加了开销。在循环外分配它们一次,然后更新它们的值。
最后一项更改应该是循环的顺序。 Numpy 数组默认是行优先(而 MATLAB 是列优先),所以ii
作为外循环表现更好。您不会注意到您大小的数组的差异,但是嘿,为什么不呢?
总的来说,修改后的函数应该是这样的。
def t_area(x, y, surface, dx):
# I assume numx == x.shape[0]. If not, pass it as an extra argument.
tri_area = np.zeros(x.shape, dtype=float)
ap = np.zeros((3,), dtype=float)
bp = np.zeros_like(ap)
cp = np.zeros_like(ap)
dp = np.zeros_like(ap)
for ii in range(x.shape[0]-1): # do you really want range(numx-1) or just range(numx)?
for jj in range(x.shape[1]-1):
xp = x[ii,jj]
yp = y[ii,jj]
zp = surface[ii,jj]
ap[:] = (xp, yp, zp)
# get `bp`, `cp` and `dp` in a similar manner and compute `tri_area[ii,jj]`