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

Question

我有一个函数（f：黑线）在特定的小区域（导数f'：蓝线，二阶导数f''：红线).我想对这个函数进行数值积分，如果我均匀分布点（在 log-space 中），我最终会在急剧变化的区域（图中 2E15 附近）出现相当大的错误。

如何构建阵列间距，使其在二阶导数较大的区域（即采样频率与二阶导数成正比）得到很好的采样？

我恰好在使用python，但我对通用算法很感兴趣。

编辑：
1)要是还能控制采样点的个数就好了（至少粗略）
2) 我考虑过构造一个形状像二阶导数的概率分布函数并从中随机抽取---但我认为这会提供较差的收敛性，一般来说，它看起来像更确定性的方法应该是可行的。

Answer 1

假设 f'' 是 NumPy array，您可以执行以下操作

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

为了仅考虑数量级波动，可以按如下方式调整...

deltas = 1/(-np.log10(1/f''))

Answer 2

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

您的数据在对数-对数图上看起来（大致）呈线性（至少，每个部分似乎是......所以，我可能会考虑在 log-space 中进行某种整合。

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

现在，对于你的每个点，你都可以获得对数对数的斜率（和截距）space:

rise = np.diff(log_y)
run = np.diff(log_x)
slopes = rise / run

同样，可以计算截距：

# y = mx + b
# :. b = y - mx
intercepts = y_log[:-1] - slopes * x_log[:-1]

好的，现在我们在 log-log 中有一堆（直线）线 space。但是，log-log space 中的一条直线对应于实际 space 中的 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) 酷方法

目前我实施了 'adaptive refinement' 方法 inspired by hydrodynamics techniques。我有一个要采样的函数 f，我选择了一些初始采样点数组 x_i。我构建了一个"sampling"函数g，它决定了在哪里插入新的样本点。

在这种情况下，我选择 g 作为 log(f) 的斜率 --- 因为我想解决日志 中的快速变化 space。然后，我将 g 的跨度划分为 L=3 个细化级别。如果 g(x_i) 超过细化级别，则该跨度被细分为 N=2 部分，这些细分被添加到样本中并针对下一个级别进行检查。这会产生这样的结果：

灰色实心线是我要采样的函数，黑色十字线是我的初始采样点
灰色虚线线是我函数对数的导数。
彩色虚线线是我的'refinement levels'
彩色十字是我细化的采样点

这一切都显示在日志中-space。

2) 简单方法

完成 (1) 后，我意识到我可能只需要在 y 中选择最大间距，然后选择 x-spacings 即可实现。同样的，把函数平均分在y，找到对应的x个点....结果如下图：

Answer 4

一种简单的方法是将 x 轴数组拆分为三个部分，并对每个部分使用不同的间距。它将允许您保持点的总数以及绘图不同区域中所需的间距。例如：

x = np.linspace(10**13, 10**15, 100)
x = np.append(x, np.linspace(10**15, 10**16, 100))
x = np.append(x, np.linspace(10**16, 10**18, 100))

您可能想根据您的数据选择更好的间距，但您明白了。

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

Construct an array spacing proportional to a function or other array

python

algorithm

numeric

numerical-methods