Точность отношения двух чисел двойной точности - PullRequest
Новая
       59

Точность отношения двух чисел двойной точности

0 голосов
/ 08 июня 2018

Допустим, я хочу вычислить численно первую производную от f (x) = cos (x) при x = pi / 2.df / dx = -sin (x), следовательно, df / dx = -1.Для этого я использую простейшую формулу:

df / dx = (f (x + h) -f (x)) / h + O (h).

Здесь O (h)это ошибка, которая пропорциональна h и, следовательно, математически O (h) стремится к нулю, а h стремится к нулю.Я понимаю, что в компьютере все по-другому.

Должен ли я ожидать, что если я использую двойную точность, которая должна давать мне от 15 до 17 десятичных значащих цифр, я смогу приблизиться к точному результату df / dx = -1.0 к ~ 10 ^ -15?Т.е. я должен найти (df / dx) _numeric + 1.0 ~ 10 ^ -15?

Это то, что я нахожу для различных значений h:

h значений: 

[0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10]

(df / dx) _numeric: 

[-0.9999999983332231,
-0.9999999999898844,
-0.9999999999175667,
-1.0000000005838656,
-0.999999993922529,
-1.000000082740371,
-1.000000082740371]

Ожидается ли это?Зачем?Почему наилучший результат получается при h = 10 ^ -5?

1 Ответ

0 голосов
/ 08 июня 2018

Расчет производных с использованием конечных разностей склонен к потере значимости .Проблема заключается в контексте, в котором вы вычитаете, а не в делении, которое вы делаете.В принципе, cos(x)-cos(x+delta) будет нормально работать, когда cos(x) и cos(x+delta) оба почти равны нулю ... но, поскольку они удаляются от нуля (и все же близки к друг другу ), точностьрезультат резко упадет.В определенный момент увеличение точности от использования меньшей дельты будет компенсировано снижением точности, вызванным потерей значимости.Для вас, похоже, что это произошло около 1e-5, но это не принципиально (и в этом регионе ошибка будет иметь тенденцию многократно колебаться).

Много было написано о том, как сделать это численностабильное конечное различие, но правило номер один - «не становитесь слишком жадным из-за малости ваших различий».Для получения более подробной (и менее расплывчатой) информации я могу порекомендовать Численные методы Формана Актона, которые (обычно) работают .

РЕДАКТИРОВАТЬ:

Извините, я забылупомянуть самый важный бит.Потеря значимости в числителе важна, но потеря значимости в x+h тем более.И эту часть на самом деле довольно легко исправить.

По мере увеличения x аппроксимация x+h ухудшается (опять же, потеря значимости).По сути, размер шага, который вы используете для оценки, не соответствует размеру шага в знаменателе.Но!Поскольку вас не особо волнует значение , точное , равное h, вы можете просто выяснить, какое значение h вы в итоге использовали после округления x+h, и использовать его в знаменателе какЧто ж.По сути, вы рассчитываете x2=x+h, а затем h'=x2-x.Вычисление h' является точным при x >> h (теорема Стербенса), исключая эту конкретную потерю значимости.

Пример кода: 

import math

def calcDerivAt_orig(f, x, h):
    x1 = x
    x2 = x+h
    y1 = f(x1)
    y2 = f(x2)
    return (y2-y1)/h

def calcDerivAt_fixed(f, x, h):
    x1 = x
    x2 = x+h
    y1 = f(x1)
    y2 = f(x2)
    return (y2-y1)/(x2-x1)

for h in [0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10]:
    dOrig = calcDerivAt_orig(math.cos, math.pi/2, h)
    origErr = abs(-1 - dOrig)
    dFixed = calcDerivAt_fixed(math.cos, math.pi/2, h)
    fixedErr = abs(-1 - dFixed)
    print("h = {}: origErr = {}, fixedErr = {}".format(h, origErr, fixedErr))

, который дает: 

h = 0.0001: origErr = 1.66677693869e-09, fixedErr = 1.66666680457e-09
h = 1e-05: origErr = 1.01155750443e-11, fixedErr = 1.6666779068e-11
h = 1e-06: origErr = 8.24332824223e-11, fixedErr = 1.66644475996e-13
h = 1e-07: origErr = 5.83865622517e-10, fixedErr = 1.55431223448e-15
h = 1e-08: origErr = 6.07747097092e-09, fixedErr = 0.0
h = 1e-09: origErr = 8.27403709991e-08, fixedErr = 0.0
h = 1e-10: origErr = 8.27403709991e-08, fixedErr = 0.0

Совсем неплохо.Конечно, мы обманываем с нашим выбором x=pi/2, потому что cos(x) там близко к нулю;на x=1.5 или что-то еще вы увидите всплывающее сообщение об ошибке при уменьшении h из-за потери значимости, которую я первоначально описал: 

h = 0.0001: origErr = 3.53519750529e-06, fixedErr = 3.5351976152e-06
h = 1e-05: origErr = 3.53675653653e-07, fixedErr = 3.53669118991e-07
h = 1e-06: origErr = 3.52831192041e-08, fixedErr = 3.53651797846e-08
h = 1e-07: origErr = 4.03034106089e-09, fixedErr = 3.44793649187e-09
h = 1e-08: origErr = 5.5453325265e-09, fixedErr = 5.1691428915e-10
h = 1e-09: origErr = 7.21702790862e-08, fixedErr = 1.03628251535e-08
h = 1e-10: origErr = 1.6659127966e-08, fixedErr = 6.58739718329e-08
...