3

以下内容来自维基百科关于牛顿方法的文章的伪代码：

#! /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 中是否有一种自动计算某种形式的fprime给定形式的方法？f

4

3 回答 3

2

f逼近at导数的一种常用方法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 时，差值接近导数，但如果太小f，该方法会由于灾难性抵消而损失精度。hx+h 也应该与 x 不同。类似的东西h = x*1e-15可能适合您的应用程序。另请参阅在 C/C++ 中实现导数。

您可以避免使用割线法逼近 f' 。它的收敛速度不如牛顿的快，但它的计算成本更低，并且您避免了必须计算导数的问题。

于 2013-05-20T14:14:39.247 回答
1

您可以近似fprime任意数量的方法。最简单的方法之一是：

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

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

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

于 2013-05-20T13:48:21.723 回答
0

回答

在. 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))

通过插入函数的系数以增加功率的顺序重新定义f使用。formula

f = formula(-2, 0, 1)

重新定义fprime，以便使用函数derivative和formula.

fprime = formula(*derivative(f))

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

概括

这是在自动计算时产生原始答案的最终解决方案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)
于 2013-05-20T15:02:02.000 回答