как использовать scipy.integrate, чтобы получить объем усеченной сферы?

Я борюсь с использованием 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

python numpy scipy numerical-integration
person samson    schedule 02.11.2012    source источник
comment
odeint используется для решения дифференциального уравнения, а не интеграла. Я действительно не понимаю, почему вы хотите использовать это здесь. Кроме того, вы можете упростить интеграл, проинтегрировав последние два измерения аналитически, чтобы у вас остался один интеграл.   -  person David Zwicker    schedule 02.11.2012

Ответы (1)


arrow_upward
3
arrow_downward

Для усечения по углу удобно использовать сферическую систему координат. Предполагая, что определение взято из Arkansas TU для radius (r), theta (t) и phi (p) как : < img src="https://i.stack.imgur.com/EKk0K.png" alt="иллюстрация в сферических координатах" />

Затем вы можете обрезать установку ограничений: 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**2 (см. mathworld ). Здесь я обрезаю половину сферы, чтобы продемонстрировать:

  • с 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]
person Saullo G. P. Castro    schedule 02.06.2013