Учитывая f, есть ли автоматический способ вычисления fprime для метода Ньютона?

Question

Следующее было перенесено из псевдокода из статьи Википедии о методе Ньютона:

#! /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)

Вопрос

Существует ли автоматизированный способ расчета какой-либо формы fprime учитывая некоторую форму f в Python 3.x?

Answer 1

Ответ

Определите функции 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))

Это должно решить ваше требование для автоматического расчета fprime от f в Python 3.x.

Резюме

Это окончательное решение, которое дает исходный ответ при автоматическом расчете 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)

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 а также f: математически разница приближается к производной, так как h стремится к 0, но метод страдает от потери точности из-за катастрофического аннулирования, если h очень маленький. Также x+h должен отличаться от x. Что-то вроде h = x*1e-15 может подойти для вашего приложения. Смотрите также реализацию производной в C / C++.

Вы можете избежать приближения f', используя секущий метод. Он не сходится так быстро, как у Ньютона, но в вычислительном отношении он дешевле, и вам не нужно вычислять производную.

Answer 3

Вы можете приблизить fprime любое количество способов. Одним из самых простых будет что-то вроде:

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

идея здесь в том, чтобы попробовать f вокруг точки x, Область выборки (определяется dx) должно быть достаточно маленьким, чтобы изменение f над этим регионом примерно линейно. Алгоритм, который я использовал, известен как метод средней точки. Вы можете получить более точные результаты, используя полиномиальные соответствия более высокого порядка для большинства функций, но это будет более затратным для вычисления.

Конечно, вы всегда будете более точными и эффективными, если будете знать аналитическую производную.