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

1

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

Мое мышление было таким: Поместим сферу с радиусом R так, чтобы ее центр находился в начале координат, и мы поместили цилиндр с радиусом R 'так, чтобы его ось была натянута на z из (b, 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()
  • 1
    Можете ли вы сначала попробовать установить b = 0, R = Rprim = 1? В этом случае ваша шаговая функция всегда равна 1 (сфера полностью находится внутри цилиндра), поэтому результат должен получиться равным 4 \ pi / 3.
  • 0
    Да, но я хочу решение для любых b, R и R '.
Показать ещё 2 комментария
Теги:
math
scipy
calculus

4 ответа

1
Лучший ответ

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

В основном я сформулировал пределы x как функцию от z и y и y как функцию от z. Тогда i, по существу, интегрирует f (z, y, z) = 1 по пересечению, используя пределы. Я сделал это из-за увеличения скорости, что позволило мне построить объем vs 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]
Поделиться
  • 0
    Рад видеть, что вы сможете найти лучшее решение :). К сожалению, я не понимаю, как вы могли бы обобщить этот подход на «более сложные функции» (какими бы они ни были). В любом случае «аналитические» решения превосходят, когда они вычисляются действительно надежным способом. Но они редко бывают :( FWIW, стохастические методы могут быть точно настроены далеко за пределы моей простой демонстрации (как по точности, так и по производительности выполнения), но IMHO считает, что такие методы имеют тенденцию быть очень надежными при обработке данных «реального мира» . Спасибо
  • 0
    Я уже сделал. Цилиндр представляет собой сферу / ион, попадающую в ядерное столкновение. Я использовал реализацию, как указано выше, чтобы рассчитать объем, непосредственно затронутый столкновением. Далее мне нужно было вычислить интеграл от потока входящего ядра, представленного цилиндром. Для этого я заменил delta_x на интеграл от функции flux по x. и определил поток как sqrt (R 2-x 2-z ** 2). После некоторой нормализации я рассчитал то, что хотел, не допуская ошибок. Я редко делаю такие вещи, это часть курса, который я беру. Благодарю.
3

Предполагая, что вы можете переводить и масштабировать свои данные таким образом, чтобы происхождение сферы находилось в [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
Поделиться
  • 1
    +1. Этот. Быстрая и грязная интеграция по методу Монте-Карло намного быстрее в настройке и гораздо менее подвержена ошибкам (и, вероятно, более эффективна), чем попытка интегрировать трехмерную функцию без регулярности путем дискретизации. И, прежде всего, он масштабируется лучше: точность равна sqrt (N), для того, чего вы хотите достичь, скорее всего, N ^ (1/3). Улучшение может заключаться в использовании квазислучайных чисел (например, последовательностей Соболя).
  • 0
    Интеграция MC, очевидно. мой мозг должен быть выключен сегодня. Благодарю.
Показать ещё 1 комментарий
2

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

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

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

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

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

В итоге у вас будет точный результат, при этом почти не требуется время вычисления.

Поделиться
  • 0
    IIRC существует умная замена переменной для пересечения шара (не сфера между прочим!) И цилиндра.
  • 0
    Оппс, я перечитал ваш комментарий и исказил его в моих предыдущих комментариях, извините.
0

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

Я бы сгенерировал тетраэдрическую сетку для домена, а затем объединил объемы ячеек. Пример с 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))

Это дает вам

Изображение 174551

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

Поделиться

Ещё вопросы

Сообщество Overcoder
Наверх
Меню