f=1/r 在 3D 任意三角形上的积分

Question

问题：

我想对函数 f=1/||x-y|| 进行积分在欧几里德 3D 空间中由一些顶点 v1,v2,v3 定义的三角形上。奇点位于 x。

问题：

如何在 python 中有效地计算任意三角形和奇点位置的积分？ scipy合适吗？

import numpy as np

# position of the singularity
x = np.array([0, 0, 1])

# triangle vertices
v1 = np.array([0, 1, 0])
v2 = np.array([1, 0, 0])
v3 = np.array([0, 0, 0])

# example for the evaluation of the function for the point pair x, y
y = np.array([0, 1, 0])
val = 1/np.linalg.norm(x-y)

# compute integral on triangle
integral = ?

Answer 1

使用纯numpy

import numpy as np

# position of the singularity
x = np.array([0, 0, 1])

# triangle vertices
v1 = np.array([0, 1, 0])
v2 = np.array([1, 0, 0])
v3 = np.array([0, 0, 0])

# compute integral on triangle
# Number of discretisation point on each axis (there are 2 axis for the surface: v2-v1 and v3-v1)
ns=1000
# I represent coordinates in triangle system: a point u,v is v1+u(v2-v1)+v(v3-v1)
# So, v1 is (0,0). v2 is (1,0), v3 is (0,1)...
# Any pair (u,v) within [0,1]×[0,1] is inside the triangle
# To let numpy do the for loop for me (never ever pure python for loop for this kind of thing)
# I use 2 dimensions. One for u, One for v. 
# Later I'll add a 3rd for the x,y,z values
# u alone use only the first dimension
u=np.linspace(0,1,ns,0)[:,None]
# Likewise, v alone use the second
# Also, note the 4th parameter of linspace: 0 (or False)
# That mean that 1 is excluded. We don't want to count both boundaries
# infinitesimal surface starting just before 1, and ending at 1 must
# be the last one
v=np.linspace(0,1,ns,0)[None,:]

# Note that u+v>=1 we are outside the triangle
inside = (u+v)<1

# The two direction vectors of my coordinate system
v21=v2-v1
v31=v3-v1

# y is a 3d array (obtained by broadcasting)
# such as, for a scalars u,v, y[u,v] is a 1D array representing `y`
# euclidean coordinates
y = v1[None,None,:]+u[...,None]*v21[None,None,:]+v[...,None]*v31[None,None,:]
# dS is just a scalar: it is the area of an elementary paralellogram
dS = np.linalg.norm(np.cross(v21,v31))/ns**2

# So, this is the 2d array of values to integrate. One for each pair (u,v)
# (0 when (u,v) is outside the triangle)
val = inside/np.linalg.norm(x-y, axis=2)
print(val.sum()*dS)

解释在注释中（最难的部分是理解注释。除此之外，这是一个非常简单的代码。它只是对每个可能的（当然，它们的子样本。离散化）位置求和所有可能的

val

三角形。 因为我使用 v2-v1 和 v3-v1 作为坐标系，所以 dS 不仅仅是 1/ns²，而只是一个标量因子。

唯一真正的计算在最后两行（我可以将它们合并为 1 行，但我想保留你的

val

）。 它为所有 y（对于所有 (u,v)）计算

1/np.linalg.norm(x-y)

（但是它是矢量化的方式：y是一个2d点数组 - 所以是一个3d数组 - 我们让

linalg.norm

计算每个点的范数 -

axis=2

- 并为我们循环执行该操作）。 然后，同样，

for

在此类

sum

的数组上完成矢量化。请注意，这可以使用其他库进行优化。或者使用 numba 显式迭代。因为我迭代了两次太多的值：一半的 (u,v) 对，即

val

数组位于三角形之外（但对总和没有贡献，因为

y

的一半是 0）。 但是与我们自己使用 for 循环相比，这个成本实际上是微不足道的（如果我们自己编写 for 循环，我们可以有

inside

，然后是

for u in linspace(0,1)

，并且永远不会在三角形之外进行迭代。但是这样做的成本在纯Python中，集成整个平行四边形的成本要高得多（数百倍），而在numba中，这将是另一个故事）。 for v in linspace(0,1-u)

确实有这样的功能。但现在我们了解了它的作用，以及为什么这些

scipy

gfun

 确实有用，就更容易理解它了。

hfun

请注意，您仍然需要找到一种与双循环集成的方法。这就是需要选择另一个坐标系（在三角形平面中）的原因，否则你不知道如何迭代 3 个坐标以在平面上积分。如果你将我的保持在

import scipy
def func(u,v):
    y=u*v21+v*v31
    return ns*ns*dS/np.linalg.norm(x-y)

print(scipy.integrate.dblquad(func, 0, 1, 0, lambda u: 1-u))

（这对于轻松计算边界很有用），那么你需要考虑变形。

dblquad

 无法知道这对 

dblquad

 代表非方形网格中的点。因此 

(u,v)

 因子 (

ns*ns*dS

 是为了撤销 dS 计算中的 

ns*ns

。

/ns**2

 中的点是积分表面（整个平行四边形）的面积与 (u,v ) 域，因此 [0,1]×[0,1])。在你的例子中，这个比率是 1 

请注意，它的速度快得令人怀疑。由于我们传递给它的是 python 的 lambda（或函数），它应该比 numpy 版本慢。所以解释是它没有那么多迭代（编辑：我刚刚算过；它做了大约 400 次采样，而我在 numpy 中采样了 1000000 次）
结果相似但不相同。
我很确定 

ns*ns*dS

比我更聪明，并且与我的常规方法相比，找到了一些非常复杂的离散方法。 但 441 个确实不多，尤其是当要集成的功能对其不透明时。所以，我认为我的 numpy 结果更接近现实