优化 Mathematica 中的插值

Question

作为我工作的一部分，我经常需要可视化复杂的 3 维密度。我使用的一个程序套件将密度的径向分量输出为对数网格上的一组 781 个点，ri = (Rmax/Rstep)^((i-1)/(pts-1)，乘以球面谐波。对于低对称性系统，球面谐波的数量可以相当大，以确保精度，例如，一个系统需要 49 个谐波，对应于 lmax = 6。在使用 v.6 并使用 Interpolation 和设置 r 构建插值径向函数时，我将得到最多 49 个插值函数的总和，每个插值函数乘以不同的球谐函数。 = Sqrt(x^2 + y^2 + z^2)，我会在一个多小时后停止 ContourPlot3D 而没有显示任何内容，这包括减少 InterpolationOrder< /code> 和 MaxRecursion 为 1。

出现了几种替代方案：

1. 评估固定网格上的密度函数，并使用 ListContourPlot 代替。
2. 或者，对径向函数进行线性样条绘制，并使用分段将它们缝合在一起。 （这本身就存在，因为我可以使用简化来帮助降低结果函数的复杂性。）

我最终使用了两者，因为 InterpolatingFunction 在其评估中给出了明显的延迟，并且最多有 49 个插值函数来评估，任何延迟都会变得明显。此外，ContourPlot3D 使用样条线速度更快，但它并没有给我带来我想要的速度。

我坦白承认，我没有在 v.7 上尝试过插值，也没有在升级后的硬件（G4 v. Intel Core i5）上尝试过。然而，我正在寻找当前方案的替代方案；最好是我可以直接使用 ContourPlot3D 的地方。我可以尝试其他形式的样条线，例如 B-spline，并可能结合使用 UnitBox 而不是使用 Piecewise。

编辑：澄清一下，我当前的实现涉及为每个径向部分创建一阶样条线，将每个样条线乘以各自的球谐函数，求和并简化上的方程每个径向间隔，然后使用分段将它们绑定到一个函数中。因此，我的实现是半解析的，因为球谐函数是精确的，并且只有径向部分是数值的。这是我希望能够使用 ContourPlot3D 的部分原因，这样我就可以利用数据的半分析性质。值得注意的是，径向网格足够精细，可以生成径向部分的良好表示，并且可以平滑插值。虽然这给了我显着的加速，但当我编写代码时，对于我当时使用的硬件来说，速度仍然很慢。

因此，我不会使用 ContourPlot3D，而是首先生成函数（如上所述），然后在 80³ 笛卡尔网格上对其进行评估。我在 ListContourPlot3D 中使用的就是这一步的数据。由于这不是自适应网格，因此在某些地方这太当然了，而且我缺少功能。

Answer 1

如果你可以不用 Mathematica，我建议你看看 Paraview （美国政府资助的 FOSS，所有平台）我发现在可视化大量数据时它优于一切。
该软件的核心是“可视化工具包”VTK，如果需要，您可以找到/编写其他前端。

VTK/Paraview 可以处理几乎任何数据类型：结构化网格上的标量和矢量或随机点、多边形、时间序列数据等。在 Mathematica 中，我经常将网格数据转储到 VTK 旧格式 在最简单的情况下看起来像这样的

# vtk DataFile Version 2.0
Generated by mma via vtkGridDump

ASCII

DATASET STRUCTURED_POINTS
DIMENSIONS 49 25 15
SPACING 0.125 0.125 0.0625
ORIGIN 8.5 5. 0.7124999999999999

POINT_DATA 18375
SCALARS  RF_pondpot_1V1MHz1amu  double 1
LOOKUP_TABLE default

0.04709501616121583
0.04135197485227461
... <18373 more numbers> ...

HTH！

Answer 2

如果确实是径向函数的插值拖慢了您的速度，您可以考虑根据您对样本点的了解对该部分进行手动编码。如下所示，这提供了显着的加速：

我用您的符号进行了设置。 lookuprvals 是一个包含 100000 个用于查找计时的 r 值的列表。

首先，以股票插值为基准

With[{interp=Interpolation[N@Transpose@{rvals,yvals}]},
  Timing[interp[lookuprvals]][[1]]]
Out[259]= 2.28466

切换到 0 阶插值已经快了一个数量级（一阶几乎相同的速度）：

With[{interp=Interpolation[N@Transpose@{rvals,yvals},InterpolationOrder->0]},
  Timing[interp[lookuprvals]][[1]]]
Out[271]= 0.146486

我们可以通过直接计算指数来获得另一个 1.5 个数量级：

Module[{avg=MovingAverage[yvals,2],idxfact=N[(pts-1) /Log[Rmax/Rstep]]},
  Timing[res=Part[avg,Ceiling[idxfact Log[lookuprvals]]]][[1]]]
Out[272]= 0.006067

作为中间立场，做手动对数线性插值。这比上面的解决方案慢，但仍然比库存插值快得多：

Module[{diffs=Differences[yvals],
  idxfact=N[(pts-1) /Log[Rmax/Rstep]]},
  Timing[Block[{idxraw,idxfloor,idxrel},
    idxraw=1+idxfact Log[lookuprvals];
    idxfloor=Floor[idxraw];
    idxrel=idxraw-idxfloor;
    res=Part[yvals,idxfloor]+Part[diffs,idxfloor]idxrel  
  ]][[1]]]
Out[276]= 0.026557

如果您有足够的内存，我会在整个网格上缓存球谐函数和半径（甚至半径索引）。然后展平网格缓存，以便您可以按

 Sum[ interpolate[yvals[lm],gridrvals] gridylmvals[lm], {lm,lmvals} ]

此处所述重新创建网格。