构造与函数或其他数组成比例的数组间距

Question

我有一个函数（f：黑线），它在一个特定的小区域内急剧变化（导数f'：蓝线，二阶导数f''：红线）。我想对这个函数进行数值积分，如果我均匀分布点（在对数空间中），那么在急剧变化的区域（图中靠近2E15）就会有相当大的误差。在

如何构造一个阵列间距，使其在二阶导数较大的区域（即与二阶导数成比例的采样频率）进行非常好的采样？

我碰巧在用python，但我对一个通用算法感兴趣。在

---

编辑：
1）如果仍然能够控制采样点的数量（至少大致如此）。
2）我曾考虑过构造一个形状类似二阶导数的概率分布函数，并从中随机抽取——但我认为这会导致收敛性较差，而且一般来说，似乎更具确定性的方法应该是可行的。

Answer 1

假设f''是NumPyarray，可以执行以下操作

# Scale these deltas as you see fit
deltas = 1/f''
domain = deltas.cumsum()

为了只考虑数量级的波动，这可以调整如下。。。在

^{pr2}$

Answer 2

我只是在这里吐口水。。。（因为我没有时间来真正的尝试一下）。。。在

你的数据看起来（大致上）在对数图上是线性的（至少，每一段看起来都是。。。所以，我可以考虑在日志空间中进行某种集成。在

log_x = log(x)
log_y = log(y)

现在，对于每个点，您可以获得log log空间中的斜率（和截距）：

^{pr2}$

同样，可以计算截距：

^{3}$

好的，现在我们在log log log空间中有一堆（直线）线。但是，log log空间中的一条直线对应于实数空间中的y = log(intercept)*x^slope。我们可以很容易地集成：y = a/(k+1) x ^ (k+1)，所以。。。在

def _eval_log_log_integrate(a, k, x):
    return np.log(a)/(k+1) * x ** (k+1)

def log_log_integrate(a, k, x1, x2):
    return _eval_log_log_integrate(a, k, x2) - _eval_log_log_integrate(a, k, x1)

partial_integrals = []
for a, k, x_lower, x_upper in zip(intercepts, slopes, x[:-1], x[1:]):
    partial_integrals.append(log_log_integrate(a, k, x_lower, x_upper))

total_integral = sum(partial_integrals)

^{你要检查一下我的数学--我已经有一段时间没有做这种事情了：-）}

Answer 3

1）酷的方法

目前，我实现了一种“自适应优化”方法inspired by hydrodynamics techniques。我有一个要采样的函数，f，我选择一些采样点的初始数组x_i。我构造了一个“sampling”函数g，它决定在哪里插入新的采样点。在

在本例中，我选择g作为log(f)的斜率，因为我想解决日志空间中的快速变化。然后我将g的范围划分为L=3细化级别。如果g(x_i)超过了一个细化级别，则该范围被细分为N=2个片段，这些细分将被添加到示例中，并根据下一个级别进行检查。结果如下：

纯灰线是我要采样的函数，黑色十字是我的初始采样点。
灰色虚线是我函数对数的导数。
彩色虚线是我的“优化级别”
彩色十字是我改进的采样点。在

这些都显示在日志空间中。在

2）简单方法

在我完成（1）之后，我意识到我可能只需要在y中选择一个最大间距，然后选择x-间距来实现这一点。类似地，只需在y中平均划分函数，并找到相应的x点。。。。结果如下：