2016-12-23 4 views
1

У меня есть эта функция (Ньютона-Рафсона алгоритм): "цифры": желательно точность корняPython Ньютона Рафсона Десятичные Корни

from scipy.misc import derivative 

def newtonDigits(function,xstart,digits): 


    xprev=0 
    ncalls=0 

    while abs(xstart-xprev) >= 0.5 * 10**(-digits): 

     ncalls +=1 
     print xstart 

     x = xstart - (function(xstart)/derivative(function,xstart,dx=1e-6)) 


     xprev = xstart 
     xstart = x 


    return xstart,ncalls 

И вход является:

f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7 *(x**3) + 20*(x**2) - 26*x + 12 

root = newtonDigits(f,1.9,6) 
print "Root: {0:.6f}".format(root[0]) 
print "Number of loops: N=" + str(root[1]) 

И выход :

1.9 
1.9347284759 
1.9570567413 
1.97161234354 
1.98117864941 
1.98749755734 
1.99168484127 
1.99446528727 
1.99631400887 
1.99754426317 
1.9983635866 
1.99890940938 
1.9992730728 
1.99951577542 
1.99967709739 
1.99978645669 
1.99985605669 
1.99990268169 
1.9999354436 
1.9999606936 
1.9999726936 
1.9999806936 
1.9999846936 
newtonr.py:55: RuntimeWarning: divide by zero encountered in double_scalars 
    x = xstart - (function(xstart)/derivative(function,xstart,dx=1e-6)) 
mainfile.py:6: RuntimeWarning: invalid value encountered in double_scalars 
    f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7 *(x**3) + 20*(x**2) - 26*x + 12 
inf 
Root: nan 
Number of loops: N=24 

Я попытался использовать Decimal, но у меня проблема с поплавками и функцией лямбда. Любые идеи?

Заранее спасибо

ответ

2

Ньютона-Рафсона не работает для корня, где производная равна нулю. Это относится к вашей проблеме, как показывает этот график.

enter image description here

В самом деле, любой корень, где производная равна нулю, трудно получить с большой точностью, и что имеет место для любого числового метода корневого ознакомительной. Ваш случай еще хуже, так как первая и вторая производные равны нулю в корне. Это означает, что график функции очень близок к прямой горизонтальной линии в корне. Учитывая это, будет трудно получить лучшую точность, чем треть общего количества значимых цифр. Большинство компьютеров будут давать двойную точность, около 15 цифр, так что более 5 значащих цифр будут сложными. Ваша процедура уже дает это, поэтому не ожидайте намного лучшего от любой цифровой процедуры.

Вы особенно беспокоитесь о том, что точность производной в вашей производной до 1e-6 (параметр в вашем вызове derivative()), и проблемы возникают, когда ваше значение x достигает одинаковой точности по x.

Вы можете использовать некоторый вариант Ньютона-Рафсона, чтобы избежать этой и другой проблемы - помните, что Ньютон-Рафсон не может сходиться в общем случае и имеет множество трудностей. Взгляните на rtsafe в книге Численные рецепты как безопасное использование N-R (Ньютон-Рафсон).

Если вы хотите придерживаться очень близко к N-R, вы можете изменить свой код, чтобы проверить значение производной, прежде чем вычислять следующее значение x. Просто прекратите цикл, когда производная станет равной нулю.

Если вам действительно нужна более высокая точность в случае, подобном вашему, и используйте прямой N-R, вам нужен лучший расчет производной. Вы используете числовое производное, которое, по-видимому, не очень хорошо. Эта числовая производная, вероятно, дает результат второго порядка, используя симметричную разность двух точек вокруг точки оценки, и это здесь недостаточно. Вы можете использовать более сложную производную процедуру, которая использует больше очков и дает результат более высокого порядка.Но в вашем случае функция достаточно проста, что вы можете написать свою собственную функцию, которая дает хорошую производную для вашей исходной функции. Просто используйте базовые правила исчисления до правила продукта. Фактически, любой реальной процедуре Ньютона-Рафсона должна быть задана функция, которая вычисляет производную, и вы обычно должны не использовать, если это необходимо, с использованием аппроксимированной числовой производной. Можно утверждать, что использование числовой производной означает, что вы на самом деле не делаете Newton-Raphson, и вы должны использовать другой метод для поиска корня - книга Числовые рецепты действительно делают это.

Модификация кода для использования определенной производной функции, а также очистки часть вашего стиля, дает

import math 

def newtonDigits(function, dfunction, xstart, digits): 
    xprev=0 
    ncalls=0 
    while abs(xstart-xprev) >= 0.5 * 10**(-digits): 
     ncalls +=1 
     print(xstart) 
     x = xstart - function(xstart)/dfunction(xstart) 
     xprev = xstart 
     xstart = x 
    return xstart, ncalls 

f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7*x**3 + 20*x**2 - 26*x + 12 
df = lambda x: 14*x*(math.e)**(x-2) + 2*(math.e)**(x-2) - 21*x**2 + 40*x - 26 

root = newtonDigits(f, df, 1.9, 6) 
print("Root: {0:.6f}".format(root[0])) 
print("Number of loops: N=" + str(root[1])) 

который печатает

1.9 
1.9347284749567717 
1.9570567399626235 
1.9716123428786936 
1.9811786535860885 
1.987497587688875 
1.991684850749086 
1.9944652840851569 
1.9963140403238206 
1.9975443983678638 
1.9983636880398847 
1.998909460162945 
1.9992731216554642 
1.9995154801019759 
1.9996770218573676 
1.999784687490283 
1.9998564592822783 
1.9999043185996854 
1.9999363386081155 
1.9999575985093117 
1.9999719257797395 
1.999982068293325 
1.9999875928805002 
2.00000490230375 
Root: 2.000005 
Number of loops: N=24 
+0

Да, проблемы возникают, когда значение x достигает такой же точности по x. Но я считаю, что это не то же самое x следующего x .. Я хочу больше точности. – tasosxak

+0

Я изменил «1e-6» на «1e-5», а результат - 1.99999953083 – tasosxak

+0

. Чистая производная подпрограмма является изворотливой и может дать странная точность, благодаря сочетанию ошибки аппроксимации, которая сжимается при более низкой допускности и ошибке округления, которая растет с меньшим допуском. Смотрите мое обширное дополнение для лучшего способа улучшить вашу рутину. Однако имейте в виду, что любой корень, где производная равна нулю, трудно получить с большой точностью, и это выполняется для любого метода поиска корней. Моя точка зрения об избежании численного дифференцирования для Ньютона-Рафсона держится. –

0

Если вы печатаете значение x и derivative на каждом шагу вы увидите:

1.9347284759 0.0683501539811 
1.95705673873 0.030810070939 
1.97161234728 0.0138147377982 
... 
1.99997209174 1.42108547152e-08 
1.99997909174 7.1054273576e-09 
inf 0.0 
nan nan 

Мы можем видеть, что производная дает нуль, и деление на ноль дает мы инф.

Если вы используете сходимость большего размера dx. Для dx=1e-3

... 
1.99991975788 2.71786149142e-06 
1.99992026519 2.7172077921e-06 
1.99992075954 2.7165576455e-06 
Root: 1.999921 
Number of loops: N=95 
+0

Ok, но почему производная равна нулю? я считаю, что это не ноль, может быть, точность не очень хорошая ... – tasosxak

Смежные вопросы