Как интегрировать 3D-данные с нерегулярной выборкой по одной оси в Python ⇐ Python

[Anonymous]

У меня есть внешняя программа (ComSol), которая предоставляет мне список образцов некоторой функции f(x,y,z).
Я хотите вычислить интеграл от f(x,y,z) по оси z, $\int_ a^b f(x,y,z) dz$.
Предоставленные выходные данные структурированы как четыре вектора: x,y,z, которые соответствуют точкам выборки, и вектор vals, который является значением функции в этих точках.
Выборки неоднородны, плотность выборки увеличивается в местах, где функция быстро меняется, и снижается там, где функция приблизительно постоянна.
Следующий код создает эквивалентный сценарий с некоторыми f(x,y,z) как некоторая базовая функция:

Код: Выделить всё

import numpy as np

def f(x,y,z) :
"function to be integrated"
return x**2 + y**2 + z**2

#generate some random non-uniform sampling points
n_points = 10000
x = np.random.random_sample(n_points)
y = np.random.random_sample(n_points)
z = np.random.random_sample(n_points)

vals = f(x,y,z)

Моя попытка решить эту проблему заключалась в интерполяции предоставленных данных в однородную сетку, которую я затем мог просто интегрировать, используя метод Симпсона.

Код: Выделить всё

import scipy
import matplotlib.pyplot as plt

#define a uniform grid to interpolate the data onto
n_interp_points = 25
x_interp = np.linspace(0,1,n_interp_points)
y_interp = np.linspace(0,1,n_interp_points)
z_interp = np.linspace(0,1,n_interp_points)

X_interp, Y_interp, Z_interp = np.meshgrid(x_interp,y_interp,z_interp)

#scipy linear interpolator with triangular basis function
tri = scipy.spatial.Delaunay(np.stack([x,y,z],axis=1)) #precompute triangles

interpolator = scipy.interpolate.LinearNDInterpolator(tri,vals,fill_value=0) #create interpolator function
vals_interp = interpolator((X_interp,Y_interp,Z_interp)) #interpolate onto defined grid

#integrate the interpolated function along the z axis with simpson's method
result = scipy.integrate.simpson(y=vals_interp,x=z_interp,axis=2)

#make comparison plots
fig, axes = plt.subplots(ncols=2,constrained_layout=True,)
for ax in axes.ravel() :
ax.set_aspect('equal')

#plot expected result
def integral_f_dz(x,y,a,b) :
"analytical expression of definite integral of f through z axis"
return 1/3 * ( -a**3 + 3*(b-a)*(x**2+y**2)+b**3 )

xx, yy = np.meshgrid(x_interp, y_interp)

axes[0].set_title("Analytically generated\nresult")
axes[0].contourf(xx, yy, integral_f_dz(xx,yy,0,1))

#plot the result from interpolation and integration
axes[1].set_title("Result from interpolation\nand integration")
axes[1].contourf(xx, yy, result)

Результирующие графики выглядят следующим образом:


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

Подробнее здесь: https://stackoverflow.com/questions/792 ... -in-python