Python - 使用 scipy 查找函数的零

Question

我很难找到（我认为）具有以下形式的相对简单函数的零（另请检查附图）：

    import numpy as np
    from numpy.random import RandomState
    from scipy import stats, optimize

    state = RandomState()

    def objective_function(x):
        initial = np.ones(10000)
        withdrawal = np.round(x, 4)

        for idx in range(175):
            sim_path = state.normal(1.001275, 0.004267, 10000)
            initial = initial - withdrawal
            initial[initial < 0] = 0
            initial = initial * sim_path

        percentile_cleared = 10000-sum(initial > 0)
        return (5000-percentile_cleared) / 10000

我一直在尝试以最少的输入使用 scipy.optimize.newton：

solution = optimize.newton(objective_function, x0=0.0075)

但实际上我很惊讶它对提供的 x0 如此敏感。 x0 的微小差异实际上决定了是否可以找到解决方案。在该特定情况下，解接近 0.0064，但我没有一个好的方法来确定它。即使在这里，提供 0.006 的 x0 也不会得到正确的解决方案。

您知道我是否应该向求解器提供更多输入，或者以不同的方式表达我的函数以使求解器更容易？提前非常感谢！

Answer 1

首先，请注意，由于

x

，您的目标函数在

np.round(x, 4)

中不可微，而牛顿方法假设该函数是可微的。然而，即使你的函数是可微的，算法也很难收敛，因为除了非常小的间隔之外，导数接近于零。

长话短说，使用 scipy.optimize.root_scalar 来查找标量函数的根。当传递包含根的区间时，它默认使用无导数方法，即

from scipy.optimize import root_scalar

root_scalar(objective_function, bracket=[0.0, 0.1])

产量

      converged: True
           flag: 'converged'
 function_calls: 38
     iterations: 36
           root: 0.006350000000384172