3d-интеграл, питон, интегрированный набор ограничений

Question

Я хотел вычислить объем пересечения сферы и бесконечного цилиндра на некотором расстоянии b, и я подумал, что сделаю это, используя быстрый и грязный скрипт на python.Мои требования - вычисление <1 с> 3 значащими цифрами.

Мое мышление было таким: мы помещаем сферу с радиусом R так, чтобы ее центр находился в начале координат, и помещаем цилиндр с радиусом R 'так, чтобы ее ось была натянута на z от(б, 0,0).Мы интегрируем по сфере, используя ступенчатую функцию, которая возвращает 1, если мы находимся внутри цилиндра, и 0, если нет, таким образом, интегрируя 1 по множеству, ограниченному нахождением внутри сферы и цилиндра, то есть пересечения.

Я попробовал это с помощью scipy.intigrate.tplquad.Это не сработало.Я думаю, что это из-за прерывания функции шага, поскольку я получаю предупреждения, такие как следующие.Конечно, я могу просто сделать это неправильно.Предполагая, что я не совершил какую-то глупую ошибку, я мог бы попытаться сформулировать диапазоны пересечения, таким образом устраняя необходимость в функции шага, но я решил, что сначала я мог бы попытаться получить некоторую обратную связь.Может ли кто-нибудь заметить какую-либо ошибку или указать на какое-то простое решение.

Предупреждение. Максимальное количество подразделений (50) достигнуто.
Если увеличение лимита не приводит к улучшению, рекомендуетсяпроанализировать
подынтегральное выражение, чтобы определить трудности.Если можно определить положение локальной сложности (особенность, разрыв), то, вероятно, выиграет от разбиения интервала и вызова интегратора на поддиапазонах.Возможно, следует использовать специальный интегратор.

Код:

from scipy.integrate import tplquad
from math import sqrt


def integrand(z, y, x):
    if Rprim >= (x - b)**2 + y**2:
        return 1.
    else:
        return 0.

def integral():
    return tplquad(integrand, -R, R, 
                   lambda x: -sqrt(R**2 - x**2),          # lower y
                   lambda x: sqrt(R**2 - x**2),           # upper y
                   lambda x,y: -sqrt(R**2 - x**2 - y**2), # lower z
                   lambda x,y: sqrt(R**2 - x**2 - y**2), # upper z
                   epsabs=1.e-01, epsrel=1.e-01
                   )

R=1
Rprim=1
b=0.5
print integral()

Answer 1

Если вы в состоянии перевести и масштабировать свои данные таким образом, что начало сферы находится в [0, 0, 0], а ее радиус - 1, тогда простое стохастическое приближение может дать вам достаточно быстрый разумный ответ.Итак, что-то в этом роде может быть хорошей отправной точкой:

import numpy as np

def in_sphere(p, r= 1.):
    return np.sqrt((p** 2).sum(0))<= r

def in_cylinder(p, c, r= 1.):
    m= np.mean(c, 1)[:, None]
    pm= p- m
    d= np.diff(c- m)
    d= d/ np.sqrt(d** 2).sum()
    pp= np.dot(np.dot(d, d.T), pm)
    return np.sqrt(((pp- pm)** 2).sum(0))<= r

def in_sac(p, c, r_c):
    return np.logical_and(in_sphere(p), in_cylinder(p, c, r_c))

if __name__ == '__main__':
    n, c= 1e6, [[0, 1], [0, 1], [0, 1]]
    p= 2* np.random.rand(3, n)- 2
    print (in_sac(p, c, 1).sum()/ n)* 2** 3

Answer 2

Выполнение тройных адаптивных числовых интегрирований над прерывной функцией, которая постоянна в двух областях, является ужасно плохой идеей, особенно если вы хотите увидеть скорость или точность.

Я бы предложил гораздо лучшую идею:чтобы уменьшить проблему аналитически.

Выровняйте цилиндр с осью путем преобразования.Это переводит сферу в некоторую точку, которая не находится в начале координат.

Теперь найдите пределы пересечения сферы с цилиндром вдоль этой оси.

Интегрируйте по этой переменной оси.Область пересечения при любом фиксированном значении вдоль оси - это просто область пересечения двух окружностей, которая, в свою очередь, просто вычисляется с использованием тригонометрии и небольших усилий.

В конце вы получите точный результат, почти не требующее времени вычислений.

Answer 3

Я решил это, используя простую интеграцию с MC, как предложено eat, но моя реализация была медленной. Мои требования возросли. Поэтому я переформулировал проблему математически, как это было предложено древесными щепками.

В основном я сформулировал пределы x как функции z и y, а y как функции z. Тогда я, по существу, интегрировал f (z, y, z) = 1 по пересечению, используя пределы. Я сделал это из-за увеличения скорости, что позволило мне построить график объема против b, и потому, что это позволило мне интегрировать более сложные функции с относительной незначительной модификацией.

Я включаю свой код на случай, если кому-то будет интересно.

from scipy.integrate import quad
from math import sqrt
from math import pi

def x_max(y,r):
    return sqrt(r**2-y**2)

def x_min(y,r):
    return max(-sqrt(r**2 - y**2), -sqrt(R**2 - y**2) + b) 

def y_max(r):
    if (R<b and b-R<r) or (R>b and b-R>r):
        return sqrt( R**2 - (R**2-r**2+b**2)**2/(4.*b**2) )
    elif r+R<b:
        return 0.
    else: #r+b<R
        return r

def z_max():
    if R>b:
        return R
    else:
        return sqrt(2.*b*R - b**2) 

def delta_x(y, r):
    return  x_max(y,r) - x_min(y,r)

def int_xy(z):
    r = sqrt(R**2 - z**2)
    return quad(delta_x, 0., y_max(r), args=(r))

def int_xyz():
    return quad(lambda z: int_xy(z)[0], 0., z_max())

R=1.
Rprim=1.
b=0.5
print 4*int_xyz()[0]

Answer 4

Первый раз: вы можете рассчитать объем пересечения вручную. Если вы не хотите (или не можете) сделать это, вот альтернатива:

Я бы сгенерировал тетраэдрическую сетку для домена и затем сложил бы объемы ячеек. Пример с pygalmesh и voropy (оба написаны мной):

import pygalmesh
import voropy

ball = pygalmesh.Ball([0, 0, 0], 1.0)
cyl = pygalmesh.Cylinder(-1, 1, 0.7, 0.1)
u = pygalmesh.Intersection([ball, cyl])

pygalmesh.generate_mesh(
        u, 'out.mesh', cell_size=0.05, edge_size=0.1
        )

mesh, _, _, _ = voropy.read('out.mesh')
print(sum(mesh.cell_volumes))

Это дает вам

и печатает 2.65692965758 как объем. Уменьшите размеры ячеек или ребер для более высокой точности.