Je me demandais si quelqu'un a su d'un package python basé sur Numpy / Scipey pour intégrer numériquement une fonction numérique compliquée sur un domaine Tessellé (dans mon cas particulier, un domaine 2D délimité par une cellule Voronoi)? Dans le passé, j'ai utilisé deux paquets de l'échange de fichiers MATLAB, mais j'aimerais rester dans mon flux de travail en Python actuel si possible. Les routines MATLAB étaient
http://www.mathworks.com/matlabentral/fileeexchange/9435-n -dimensionale-simplex-quadrature
pour la génération quadrature et maille utilisant:
http://www.mathworks.com/matlabentral/fileExchange/25555-mesh2d -Automatic-maillage-génération
Toute suggestion sur la génération de maille puis une intégration numérique sur ce maillage serait appréciée.
3 Réponses :
Je ne suis pas sûr que dblquad gère facilement le domaine de l'intégration dont j'ai besoin.
On dirait que vous avez une grille assez désordonnée et discrète que vous intégrez. La communauté des éléments finis a probablement quelque chose pour vous, mais ce n'est pas mon domaine d'expertise.
Cela intègre directement les triangles directement, pas les régions de Voronoi, mais devraient être proches. (Exécutez avec différents nombres de points à voir?) Cela fonctionne également dans 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
Merci pour votre réponse. Je vais tester cela un point bientôt pour voir si cela fonctionne pour mon cas d'utilisation.
@JOSH, que diriez-vous d'une fonction de test avec une intégrale connue, de comparer des méthodes?
L'intégration numérique est typiquement définie sur des éléments simples tels que les triangles ou les rectangles. Cela dit, vous pouvez décomposer chaque polonge dans une série de triangles. Avec une chance, votre maillage polygonal a un double triangulaire qui faciliterait beaucoup plus facilement l'intégration.
quadpy (un projet de mien) L'intégration numérique sur de nombreux domaines, par exemple des triangles: Vous pouvez également intégrer de manière non adaptative en fournissant une de centaines de schémas pour le triangle.
Quadpy est un package merveilleux pour l'intégration spécifique d'application, en particulier pour la méthode des éléments finis, merci!
Votre entrée est-elle dispersée Data XI, Yi, ZI ou est-ce (comme vous écrivez) une seule cellule Voronoi à la fois, avec ~ 6 côtés? Aussi, combien de points ou de cellules - 1000, 1000000?
@Denis: C'est une collection de n points XI, yi qui marquent les centres des cellules Voronoi, où n est de l'ordre de 10 ^ 2 à 10 ^ 3. Le nombre de côtés de chaque cellule Voronoi n'est pas garanti d'être un nombre particulier.
Josh, que diriez-vous des tags Mesh Qhull (- Python)?
Consultez Quadrature . Vous trouverez que de nombreux régimes pour triangles et Tets sont ici. Pour la génération de mailles, il y a MSHR , pygmsh , MSHPY , Frentos etc.