给定f，牛顿法有没有自动计算fprime的方法？

Question

以下内容是从维基百科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)

问题

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

Answer 1

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

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

这里的想法是围绕点f取样x。取样区域（由dx决定）应足够小，使{}在该区域上的变化近似为线性。我使用的算法称为中点法。你可以通过对大多数函数使用高阶多项式拟合来获得更精确的结果，但这将是更昂贵的计算。在

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

Answer 2

近似f在x处的导数的一种常见方法是使用有限差分：

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

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

可以使用secant method来避免近似f'。它的收敛速度没有牛顿的快，但它的计算成本更低，而且避免了计算导数的问题。在

Answer 3

回答

将函数formula和derivative直接定义为import之后。在

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，方法是按照幂递增的顺序插入函数的系数。在

^{pr2}$

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

fprime = formula(*derivative(f))

这应该可以解决Python3.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)