4

scipy.integrate の使用に苦労しています。tplquad を使用integrateしましたが、(切り捨てられた) 球体のボリュームを取得するにはどうすればよいですか? ありがとう

import scipy
from scipy.integrate import quad, dblquad, tplquad
from math import*
from numpy import *

R = 0.025235 #radius
theta0 = acos(0.023895) #the angle from the edge of truncated plane to the center of
sphere

def f_1(phi,theta,r):
    return r**2*sin(theta)*phi**0
Volume = tplquad(f_1, 0.0,R, lambda y: theta0, lambda y: pi, lambda y,z: 0.0,lambda
y,z: 2*pi)

print Volume
4

1 に答える 1

3

角度で切り捨てるには、球座標系を使用すると便利です。、およびとしてアーカンソー工科大学から取得した定義を仮定すると、次のようになります。 radius (r)theta (t)phi (p)球座標の図

次に、制限の設定を切り捨てることができますr1 r2 t1 t2 p1 p2:

import scipy
from scipy.integrate import quad, dblquad, tplquad
from numpy import *

# limits for radius
r1 = 0.
r2 = 1.

# limits for theta
t1 = 0
t2 = 2*pi

# limits for phi
p1 = 0
p2 = pi

def diff_volume(p,t,r):
    return r**2*sin(p)

volume = tplquad(diff_volume, r1, r2, lambda r:   t1, lambda r:   t2,
                                      lambda r,t: p1, lambda r,t: p2)[0]

平面で切り捨てるには、デカルト座標系を使用すると便利です(x,y,z)( x**2+y**2+z**2=R**2mathworldを参照)。ここでは、デモンストレーションのために球体の半分を切り捨てています。

  • からx1=-Rまでx2=R
  • からy1=0までy2=(R**2-x**2)**0.5
  • からz1=-(R**2-x**2-y**2)**0.5までz2=(R**2-x**2-y**2)**0.5

ラムダを使用した便利な例:

R= 2.

# limits for x
x1 = -R
x2 = R

def diff_volume(z,y,x):
    return 1.

volume = tplquad(diff_volume, x1, x2,
                 lambda x: 0., lambda x: (R**2-x**2)**0.5,
                 lambda x,y: -(R**2-x**2-y**2)**0.5,
                 lambda x,y:  (R**2-x**2-y**2)**0.5 )[0]
于 2013-06-02T16:11:43.200 に答える