scipy 退化寻根问题:得到f(x)>0的第一个值

gkl3eglg  于 2023-05-07  发布在  其他
关注(0)|答案(2)|浏览(106)

给定函数 f(x),对于小于临界值 c 的所有值 x 为零,并且对于值 x〉c 为非零。
我想用优化方法近似计算临界值 c。因为函数 f(x) 的开销很大,所以我想尽可能少地计算它。因此,针对值 x 的预定义列表计算 f(x) 是不可行的。
考虑一个像下面这样的函数,临界点 a=1.4142

def func(x):
    return max(x**2 - 2, 0) if x > 0 else 0

我会使用一些自定义的二分函数来实现这一点。然而,我想知道这个“退化”的根查找问题是否可以使用SciPy或NumPy中的现有函数来解决。我用scipy.optimize.root_scalar做了实验。但是,它似乎不支持上面的功能。

from scipy.optimize import root_scalar
root_scalar(func, bracket=[-6, 8])
# Yields an error: f(a) and f(b) must have different signs.
yyyllmsg

yyyllmsg1#

从技术上讲,该函数中x〈1.414的任何地方都是根。如果你想使用根查找器,你需要使它成为一个有效的根查找问题。我建议从函数的输出中减去一个非常小的值:

from scipy.optimize import root_scalar
import math
wrapper = lambda func: lambda *args: func(*args) - math.ulp(0)
root_scalar(wrapper(func), bracket=[-6, 8])

这是通过从函数中减去一个ULP(最后一个单位,可以添加或减去的最小值)来实现的。这使得它对所有x 1.414都是严格负< 1.414 and strictly positive for all x >的。在一台64位浮点运算的计算机上,这意味着它从每个结果中减去5*10-324,这会稍微移动根的位置,但不会太多。

kkih6yb8

kkih6yb82#

这只是为了呈现适用于所呈现的示例的“自定义”二等分方法。但是,我仍然建议使用公认的答案。

def bisect(func, x_min, x_max, x_func=None,
           tol=None, iter_max=None, **kwargs):
    y_min = func(x_min, **kwargs)
    if y_min > 0:
        msg = "Warning: no solution as y_min>0, with x_min=%f."
        print(msg % x_min)
        return x_min
    y_max = func(x_max, **kwargs)
    if y_max <= 0:
        msg = "Warning: no solution as y_max<=0, with x_max=%f."
        print(msg % x_max)
        return x_max
    if tol is None and iter_max is None:
        tol = 1e-7
    if x_func is None:
        x_func = lambda x0, x1: (x1+x0)/2
    
    from itertools import count, islice
    x_last = np.infty
    for cnt in islice(count(), iter_max):
        x_new = x_func(x_min, x_max)
        y_new = func(x_new, **kwargs)
        if y_new<=0:
            x_min = x_new
        else:
            x_max = x_new
        if (tol is not None) and abs(x_last - x_new) < tol:
            break
        x_last = x_new
    return x_min if y_new>0 else x_max

def func(x):
    return max(x**2 - 2, 0) if x > 0 else 0

ret = bisect(func, x_min=-6, x_max=8)
print(ret)

有趣的是,对于这种基本类型的函数func(x)bisect()的性能似乎优于root_scalar方法。但我预计,对于更复杂的函数,这种优势将消失。

%timeit bisect(func, x_min=-6, x_max=8, iter_max=None, tol=1e-7)
# 46.3 µs ± 1.96 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

%timeit root_scalar(wrapper(func), bracket=[-6, 8], xtol=1e-7)
# 84.6 µs ± 6.29 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

# Test that the results are equivalent:
x1 = bisect(func, x_min=-6, x_max=8, iter_max=None, tol=1e-7)
x2 = root_scalar(wrapper(func), bracket=[-6, 8], xtol=1e-7).root
assert abs(x2-x1) < 1e-7
# 2.9078300656237843e-08

相关问题