2

我想根据三维离散不均匀密度分布计算球体的质量。假设一组不同密度的 3x3x3 立方体被一个球体内接。使用 Python 总结分区质量的最快方法是什么？

我试图计算一个球体的数学方程下的体积： x^2 +y^2 +z^2 = R^2 用于使用scipy.integrate.dblquad的立方体之一的范围。但是，只有当边界小于球体的半径时，结果才有效，并且重复计算假设 50'000 个球体，每个球体有 27 个立方体会很慢。

另一方面，由于质量分布相当粗糙和离散，我认为不能使用通常的 CoM 计算方程。

4

3 回答 3

1

计时实验

你没有指定你的时间限制，所以我用一个很好的集成包做了一个小实验。

如果立方体密度是简单的函数，则在没有优化的情况下，可以在标准笔记本电脑中在 0.005 秒内评估球坐标中的每个积分。

作为参考，这是 Mathematica 中的程序： 

Clear@f;
(* Define a cuboid as density function *)
iP = IntegerPart;
f[{x_, y_, z_}, {lx_, ly_, lz_}] :=   iP[x - lx] + iP[y - ly] + iP[z - lz] /; 
   lx <= x <= lx + 3 && ly <= y <= ly + 3 && lz <= z <= lz + 3;

f[{x_, y_, z_}, {lx_, ly_, lz_}] := Break[] /; True;

Timing[Table[s = RandomReal[{0, 3}, 3]; (*sphere center random*)
   sphereRadius = Min[Union[s, 3 - s]]; (*max radius inside cuboid *)
   NIntegrate[(f[{x, y, z} - s, -s] /.  (*integrate in spherical coords *)
       {x -> r Cos@th Sin@phi, 
        y -> r Sin@th Sin@phi, 
        z -> r Cos@phi}) r^2 Sin@phi,
       {r, 0, sphereRadius}, {th, 0, 2 Pi}, {phi, 0, Pi}], 
         {10000}]][[1]]

结果是52 秒，10^4 次迭代。

所以也许你不需要优化很多......

于 2010-12-13T00:27:36.000 回答
0

我无法得到您对球体刻字的确切含义。我也没有尝试过 scipy.integrate。但是，这里有一些：

将 3x3x3 立方体设置为单位密度。然后分别对每个立方体进行积分，所以你应该在这里有体积立方体 V_ijk。现在对于每个球体，您可以通过求和得到每个球体的质量V_ijk*D_ijk，其中D_ijk是球体的密度。

它应该快得多，因为您现在不需要进行集成。

于 2010-12-10T19:15:19.983 回答
0

您可以获得立方体（或矩形棱柱）和球体之间相交体积的解析公式。这并不容易，但应该是可能的。我已经为 2D 中的任意三角形和圆做了它。基本思想是将交叉点分解为更简单的部分，例如四面体和体积球面三角形扇区，其中相对简单的体积公式是已知的。主要困难在于考虑所有可能的交叉点情况。幸运的是，这两个对象都是凸的，因此可以保证单个凸相交体积。

一种近似方法可能是简单地细分立方体，直到您的近似数值积分算法起作用；这应该还是比较快的。你知道匹克定理吗？这仅适用于 2D，但我相信有3D 概括。

于 2010-12-10T21:44:15.393 回答