我想知道是否有人知道基于 numpy/scipy 的 python 包可以在曲面细分域（在我的特定情况下，一个由 voronoi 单元限制的二维域）上对复杂的数值函数进行数值积分？过去，我使用了 matlab 文件交换中的几个包，但如果可能的话，我想留在我当前的 python 工作流程中。matlab例程是
http://www.mathworks.com/matlabcentral/fileexchange/9435-n-dimensional-simplex-quadrature
对于正交和网格生成，使用：
http://www.mathworks.com/matlabcentral/fileexchange/25555-mesh2d-automatic-mesh-generation
任何关于网格生成以及对该网格进行数值积分的建议都将不胜感激。
这直接整合了三角形，而不是 Voronoi 区域，但应该是接近的。（运行不同数量的点来查看？）它也适用于 2d、3d ...
#!/usr/bin/env python
from __future__ import division
import numpy as np
__date__ = "2011-06-15 jun denis"
#...............................................................................
def sumtriangles( xy, z, triangles ):
""" integrate scattered data, given a triangulation
zsum, areasum = sumtriangles( xy, z, triangles )
In:
xy: npt, dim data points in 2d, 3d ...
z: npt data values at the points, scalars or vectors
triangles: ntri, dim+1 indices of triangles or simplexes, as from
http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html
Out:
zsum: sum over all triangles of (area * z at midpoint).
Thus z at a point where 5 triangles meet
enters the sum 5 times, each weighted by that triangle's area / 3.
areasum: the area or volume of the convex hull of the data points.
For points over the unit square, zsum outside the hull is 0,
so zsum / areasum would compensate for that.
Or, make sure that the corners of the square or cube are in xy.
"""
# z concave or convex => under or overestimates
npt, dim = xy.shape
ntri, dim1 = triangles.shape
assert npt == len(z), "shape mismatch: xy %s z %s" % (xy.shape, z.shape)
assert dim1 == dim+1, "triangles ? %s" % triangles.shape
zsum = np.zeros( z[0].shape )
areasum = 0
dimfac = np.prod( np.arange( 1, dim+1 ))
for tri in triangles:
corners = xy[tri]
t = corners[1:] - corners[0]
if dim == 2:
area = abs( t[0,0] * t[1,1] - t[0,1] * t[1,0] ) / 2
else:
area = abs( np.linalg.det( t )) / dimfac # v slow
zsum += area * z[tri].mean(axis=0)
areasum += area
return (zsum, areasum)
#...............................................................................
if __name__ == "__main__":
import sys
from time import time
from scipy.spatial import Delaunay
npt = 500
dim = 2
seed = 1
exec( "\n".join( sys.argv[1:] )) # run this.py npt= dim= ...
np.set_printoptions( 2, threshold=100, edgeitems=5, suppress=True )
np.random.seed(seed)
points = np.random.uniform( size=(npt,dim) )
z = points # vec; zsum should be ~ constant
# z = points[:,0]
t0 = time()
tessellation = Delaunay( points )
t1 = time()
triangles = tessellation.vertices # ntri, dim+1
zsum, areasum = sumtriangles( points, z, triangles )
t2 = time()
print "%s: %.0f msec Delaunay, %.0f msec sum %d triangles: zsum %s areasum %.3g" % (
points.shape, (t1 - t0) * 1000, (t2 - t1) * 1000,
len(triangles), zsum, areasum )
# mac ppc, numpy 1.5.1 15jun:
# (500, 2): 25 msec Delaunay, 279 msec sum 983 triangles: zsum [ 0.48 0.48] areasum 0.969
# (500, 3): 111 msec Delaunay, 3135 msec sum 3046 triangles: zsum [ 0.45 0.45 0.44] areasum 0.892
数值积分通常是在三角形或矩形等简单元素上定义的。也就是说，您可以将每个 polgon 分解为一系列三角形。运气好的话，您的多边形网格有一个三角形对偶，这将使集成更加容易。
quadpy（我的一个项目）对许多域进行数值积分，例如三角形：
import numpy
import quadpy
sol, error_estimate = quadpy.t2.integrate_adaptive(
lambda x: numpy.exp(x[0]),
numpy.array(
[
[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]],
[[1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0]],
[[0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]],
]
),
1.0e-10,
)
print(sol)
3.5914091422921017
您还可以通过为三角形提供数百种方案之一来进行非自适应集成。
scipy.integrate.dblquad怎么样？它使用自适应正交规则，因此您放弃对集成网格的控制。不知道这对您的应用程序是有利还是不利。