python数值积分_区域体积的Python数值积分

For a program, I need an algorithm to very quickly compute the volume of a solid. This shape is specified by a function that, given a point P(x,y,z), returns 1 if P is a point of the solid and 0 if P is not a point of the solid.

I have tried using numpy using the following test:

import numpy

from scipy.integrate import *

def integrand(x,y,z):

if x**2. + y**2. + z**2. <=1.:

return 1.

else:

return 0.

g=lambda x: -2.

f=lambda x: 2.

q=lambda x,y: -2.

r=lambda x,y: 2.

I=tplquad(integrand,-2.,2.,g,f,q,r)

print I

but it fails giving me the following errors:

Warning (from warnings module):

File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321

warnings.warn(msg, IntegrationWarning)

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.

If increasing the limit yields no improvement it is advised to analyze

the integrand in order to determine the difficulties. If the position of a

local difficulty can be determined (singularity, discontinuity) one will

probably gain from splitting up the interval and calling the integrator

on the subranges. Perhaps a special-purpose integrator should be used.

Warning (from warnings module):

File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321

warnings.warn(msg, IntegrationWarning)

IntegrationWarning: The algorithm does not converge. Roundoff error is detected

in the extrapolation table. It is assumed that the requested tolerance

cannot be achieved, and that the returned result (if full_output = 1) is

the best which can be obtained.

Warning (from warnings module):

File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321

warnings.warn(msg, IntegrationWarning)

IntegrationWarning: The occurrence of roundoff error is detected, which prevents

the requested tolerance from being achieved. The error may be

underestimated.

Warning (from warnings module):

File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 321

warnings.warn(msg, IntegrationWarning)

IntegrationWarning: The integral is probably divergent, or slowly convergent.

So, naturally, I looked for "special-purpose integrators", but could not find any that would do what I needed.

Then, I tried writing my own integration using the Monte Carlo method and tested it with the same shape:

import random

# Monte Carlo Method

def get_volume(f,(x0,x1),(y0,y1),(z0,z1),prec=0.001,init_sample=5000):

xr=(x0,x1)

yr=(y0,y1)

zr=(z0,z1)

vdomain=(x1-x0)*(y1-y0)*(z1-z0)

def rand((p0,p1)):

return p0+random.random()*(p1-p0)

vol=0.

points=0.

s=0. # sum part of variance of f

err=0.

percent=0

while err>prec or points

p=(rand(xr),rand(yr),rand(zr))

rpoint=f(p)

vol+=rpoint

points+=1

s+=(rpoint-vol/points)**2

if points>1:

err=vdomain*(((1./(points-1.))*s)**0.5)/(points**0.5)

if err>0:

if int(100.*prec/err)>=percent+1:

percent=int(100.*prec/err)

print percent,'% complete\n error:',err

print int(points),'points used.'

return vdomain*vol/points

f=lambda (x,y,z): ((x**2)+(y**2)<=4.) and ((z**2)<=9.) and ((x**2)+(y**2)>=0.25)

print get_volume(f,(-2.,2.),(-2.,2.),(-2.,2.))

but this works too slowly. For this program I will be using this numerical integration about 100 times or so, and I will also be doing it on larger shapes, which will take minutes if not an hour or two at the rate it goes now, not to mention that I want a better precision than 2 decimal places.

I have tried implementing a MISER Monte Carlo method, but was having some difficulties and I'm still unsure how much faster it would be.

So, I am asking if there are any libraries that can do what I am asking, or if there are any better algorithms which work several times faster (for the same accuracy). Any suggestions are welcome, as I've been working on this for quite a while now.

EDIT:

If I cannot get this working in Python, I am open to switching to any other language that is both compilable and has relatively easy GUI functionality. Any suggestions are welcome.

解决方案

As the others have already noted, finding the volume of domain that's given by a Boolean function is hard. You could used pygalmesh (a small project of mine) which sits on top of CGAL and gives you back a tetrahedral mesh. This

import numpy

import pygalmesh

import meshplex

class Custom(pygalmesh.DomainBase):

def __init__(self):

super(Custom, self).__init__()

return

def eval(self, x):

return (x[0]**2 + x[1]**2 + x[2]**2) - 1.0

def get_bounding_sphere_squared_radius(self):

return 2.0

mesh = pygalmesh.generate_mesh(Custom(), cell_size=1.0e-1)

gives you

From there on out, you can use a variety of packages for extracting the volume. One possibility: meshplex (another one out of my zoo):

import meshplex

mp = meshplex.MeshTetra(mesh.points, mesh.cells["tetra"])

print(numpy.sum(mp.cell_volumes))

gives

4.161777

which is close enough to the true value 4/3 pi = 4.18879020478.... If want more precision, decrease the cell_size in the mesh generation above.