给定 f，是否有一种自动计算牛顿法的 fprime 的方法？答案

【问题标题】：Given f, is there an automatic way to calculate fprime for Newton's method?给定 f，是否有一种自动计算牛顿法的 fprime 的方法？
【发布时间】：2013-05-15 01:54:36
【问题描述】：

以下内容来自维基百科文章Newton's method的伪代码：

#! /usr/bin/env python3
# https://en.wikipedia.org/wiki/Newton's_method
import sys

x0 = 1
f = lambda x: x ** 2 - 2
fprime = lambda x: 2 * x
tolerance = 1e-10
epsilon = sys.float_info.epsilon
maxIterations = 20

for i in range(maxIterations):
    denominator = fprime(x0)
    if abs(denominator) < epsilon:
        print('WARNING: Denominator is too small')
        break
    newtonX = x0 - f(x0) / denominator
    if abs(newtonX - x0) < tolerance:
        print('The root is', newtonX)
        break
    x0 = newtonX
else:
    print('WARNING: Not able to find solution within the desired tolerance of', tolerance)
    print('The last computed approximate root was', newtonX)

问题

在 Python 3.x 中给定某种形式的 f，是否有一种自动计算某种形式的 fprime 的方法？

【问题讨论】：

不直接回答，但最好将其保存为列表，例如 x**2 -2 = [1,0,-2] 2x 将是 [0,2,0]
对于f的一些种形式，您可以轻松自动计算出f'的对应形式。例如，如果您将多项式表示为系数列表（或映射exponent -> coefficient），则可以轻松计算导数。您可以超越这一点并包括其他基本功能，但一般情况下的答案是否定的。但是您可以使用泰勒展开式和仅f 的值来近似任何给定点的导数。
感谢 Beasley 和 Fischer 提供帮助的 cmets。另一个问题被证明是有用的：stackoverflow.com/questions/11155367

标签： python math python-3.x solver newtons-method

【解决方案1】：

您可以通过多种方式近似fprime。最简单的方法之一是：

lambda fprime x,dx=0.1: (f(x+dx) - f(x-dx))/(2*dx)

这里的想法是围绕点x 对f 进行采样。采样区域（由dx 确定）应该足够小，以使f 在该区域上的变化近似为线性。我使用的算法称为中点法。您可以通过对大多数函数使用高阶多项式拟合来获得更准确的结果，但计算成本会更高。

当然，如果你知道解析导数，你总是会更加准确和高效。

【讨论】：

也许我想要维基百科文章中伪代码的 Python 版本，然后认为如果只给定 f 变量，自动生成 fprime 变量会很好......

【解决方案2】：

在x 处逼近f 的导数的常用方法是使用有限差分：

f'(x) = (f(x+h) - f(x))/h                   Forward difference
f'(x) = (f(x+h) - f(x-h))/2h                Symmetric

h 的最佳选择取决于x 和f：从数学上讲，当 h 趋于 0 时，差值逼近导数，但如果 h 太高，该方法会由于灾难性抵消而损失精度小的。 x+h 也应该与 x 不同。 h = x*1e-15 之类的东西可能适合您的应用程序。另见implementing the derivative in C/C++。

您可以通过使用secant method 来避免逼近 f'。它的收敛速度不如牛顿的快，但计算成本更低，而且您避免了必须计算导数的问题。

【讨论】：

如果f 的形式合适，为什么要近似导数？
牛顿的方法不限于任何特定形式的函数。您想将您的实现限制为特定形式的功能吗？哪一个？（多项式，有理数，三角函数，...）
啊，我们对“形式”这个词有不同的含义。您的意思是方程的形式，而我的意思是程序表示的形式。由于给定方程是具有单个变量的多项式，因此这是我求解fprime 的目标。数学家可能会针对更通用的解决方案。
给定一个编程表示，您可以使用符号方法获得导数，例如使用 sympy：docs.sympy.org/0.7.2/tutorial.html#calculus 但是，导数的符号表示可能比近似值更慢。或者，也许一个 sympy 表达式可以编译为字节码，这样它就可以像“本机”python 表达式一样快，我不知道。

【解决方案3】：

回答

在import 之后直接定义函数formula 和derivative，如下所示。

def formula(*array):
    calculate = lambda x: sum(c * x ** p for p, c in enumerate(array))
    calculate.coefficients = array
    return calculate

def derivative(function):
    return (p * c for p, c in enumerate(function.coefficients[1:], 1))

使用formula 重新定义f，方法是插入函数的系数以增加功率。

f = formula(-2, 0, 1)

重新定义fprime，以便使用derivative 和formula 函数自动创建。

fprime = formula(*derivative(f))

这应该可以解决您在 Python 3.x 中从 f 自动计算 fprime 的要求。

总结

这是在自动计算fprime时产生原始答案的最终解决方案。

#! /usr/bin/env python3
# https://en.wikipedia.org/wiki/Newton's_method
import sys

def formula(*array):
    calculate = lambda x: sum(c * x ** p for p, c in enumerate(array))
    calculate.coefficients = array
    return calculate

def derivative(function):
    return (p * c for p, c in enumerate(function.coefficients[1:], 1))

x0 = 1
f = formula(-2, 0, 1)
fprime = formula(*derivative(f))
tolerance = 1e-10
epsilon = sys.float_info.epsilon
maxIterations = 20

for i in range(maxIterations):
    denominator = fprime(x0)
    if abs(denominator) < epsilon:
        print('WARNING: Denominator is too small')
        break
    newtonX = x0 - f(x0) / denominator
    if abs(newtonX - x0) < tolerance:
        print('The root is', newtonX)
        break
    x0 = newtonX
else:
    print('WARNING: Not able to find solution within the desired tolerance of', tolerance)
    print('The last computed approximate root was', newtonX)

【讨论】：