寻找单变量大量函数的根

Question

寻找单变量大量函数的根

pythonnumpyscipyraytracing

11

我正在使用Python/numpy/scipy编写一个小型的光线跟踪器。表面被建模为二维函数，其给出相对于法平面的高度。我将射线与表面的交点问题简化为找到一个单变量函数的根。这些函数是连续可微的。

有没有比简单遍历所有函数更有效的方法，可以使用scipy根查找器（可能还要使用多个进程）？

编辑：这些函数是表示射线和表面函数之间差异的函数，受到交点平面的限制。

- mikebravo

2

什么是函数？它可能有解析解吗？ - mgilson

1

表面函数可以任意选择，我希望它具有灵活性。对于特定的函数（即Chebyshev多项式的叠加），存在解析解，但可能涉及许多参数。通过解决线性方程组来找到交点应该对特定的表面是可行的。 - mikebravo

有标准的方法可以找到射线/平面、射线/球体、射线/三角形的交点。你能将你的表面建模成一个三角网格吗？如果没有解析解或几何逼近你的表面函数，我不知道是否有比直接计算函数更有效的方法。 - engineerC

我考虑过将函数离散化，但是表面正在被优化并不断变化 - 我希望有一个三角剖分，可以将整个区域均匀地分成面。我已经在二维等距分割方面遇到了一些困难。 - mikebravo

2个回答

1

在过去的几年中，scipy.optimize.newton 获得了矢量化支持。使用另一个答案中的示例现在看起来像这样：

import numpy as np
from scipy import optimize

def F(x, a, b):
    return np.power(x, a+1.0) - b

N = 1000000

a = np.random.rand(N)
b = np.random.rand(N)

optimize.newton(F, np.zeros(N), args=(a, b))

这个方法的运行速度和其他答案中向量化二分法一样快。

- Akira Kyle

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Alex I · Accepted Answer

以下示例展示了如何使用二分法并行计算函数 x**(a+1) - b （其中 a 和 b 不同）的一百万个根。这将花费约 12 秒钟的时间。

import numpy

def F(x, a, b):
    return numpy.power(x, a+1.0) - b

N = 1000000

a = numpy.random.rand(N)
b = numpy.random.rand(N)

x0 = numpy.zeros(N)
x1 = numpy.ones(N) * 1000.0

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F0 = F(x0, a, b)
    F1 = F(x1, a, b)
    F_mid = F(x_mid, a, b)
    x0 = numpy.where( numpy.sign(F_mid) == numpy.sign(F0), x_mid, x0 )
    x1 = numpy.where( numpy.sign(F_mid) == numpy.sign(F1), x_mid, x1 )
    error_max = numpy.amax(numpy.abs(x1 - x0))
    print "step=%d error max=%f" % (step, error_max)
    if error_max < 1e-6: break

基本思想是并行运行向量变量的所有常规根查找步骤，使用可在一组变量向量和定义各个组件函数的等效向量参数上评估函数。条件句被替换为掩码和numpy.where()的组合。这可以持续进行，直到找到所需精度的所有根，或者直到找到足够多的根，值得将其从问题中删除，并继续使用排除这些根的较小问题。

我选择要解决的函数是任意的，但如果函数行为良好，那么会有所帮助；在本例中，家族中的所有函数都是单调的，并且具有唯一的正根。此外，对于二分法，我们需要变量的猜测值，给出不同函数符号，而这些在这里也很容易想到（x0和x1的初始值）。

以上代码使用可能是最简单的根查找器（二分法），但是相同的技术可以轻松应用于牛顿-拉弗森、Ridder等。根查找方法中条件语句越少，它就越适合用于此。然而，您必须重新实现您想要的任何算法，没有办法直接使用现有的库根查找功能。

以上代码片段是以清晰度为主要考虑，而不是速度。避免某些计算的重复，特别是每次迭代只评估一次函数而不是3次，将其加速到9秒，如下所示：

...
F0 = F(x0, a, b)
F1 = F(x1, a, b)

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F_mid = F(x_mid, a, b)
    mask0 = numpy.sign(F_mid) == numpy.sign(F0)
    mask1 = numpy.sign(F_mid) == numpy.sign(F1)
    x0 = numpy.where( mask0, x_mid, x0 )
    x1 = numpy.where( mask1, x_mid, x1 )
    F0 = numpy.where( mask0, F_mid, F0 )
    F1 = numpy.where( mask1, F_mid, F1 )
...

相比之下，使用scipy.bisect()逐个查找根需要约94秒：

for i in range(N):
    x_root = scipy.optimize.bisect(lambda x: F(x, a[i], b[i]), x0[i], x1[i], xtol=1e-6)