给定一个几何对象,为简化起见,它是一个具有一定半径的半球。它显示为一个2D矩阵,Z数据为高度。假设我沿着任意线切割对象,我想计算切割的面积。我的解决方案是使用scipys RectBivariateSpline对半球进行插值,以精确显示它。
import numpy as np
import scipy.interpolate as intp
radius = 15.
gridsize = 0.5
spectrum = np.arange(-radius,radius+gridsize,gridsize)
X,Y = np.meshgrid(spectrum,spectrum)
Z = np.where(np.sqrt(X**2+Y**2)<=radius, np.sqrt(radius**2-np.sqrt(X**2+Y**2)**2), 0)
spline = intp.RectBivariateSpline(x = X[0,:], y = Y[:,0], z = Z)
# Example coordinates of the cut
x0 = -4.78
x = -6.73
y0 = -15.
y = 15.
然而,RectBivariateSpline只提供面积积分(可以通过设置x0 = x或y0 = y快速检查)。另一方面,UnivariateSpline只接受一维数组,只有当我的切割碰巧沿着矩阵Z的一个特定向量时才能工作。
因为我想把这个运算进行几千次,所以我需要一个相对快速的方法来解这个积分(只要误差可以忽略不计,数值或解析都无所谓)。有人知道怎么做吗?
1条答案
按热度按时间rur96b6h1#
事实证明,对于我的情况,沿着我的切割对样条进行采样(使用numpy的arange来收集等间距的点),然后通过辛普森规则进行积分就足够了,辛普森规则只需要距离足够小的点(可以通过arange的step参数来控制)。