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

首先，请注意，由于np.round(x, 4)，您的目标函数在x 中不可微分，而牛顿法假定该函数是可微分的。然而，即使你的函数是可微的，算法也很难收敛，因为除了非常小的区间外，导数接近于零。

长话短说，使用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