Числовое дифференцирование Python и минимальное значение для h

Question

Я вычисляю первую производную, используя следующий код:

def f(x):
   f = np.exp(x)
   return f

def dfdx(x):
   Df = (f(x+h)-f(x-h)) / (2*h)
   return Df

Например, для x == 10 это работает нормально.Но когда я устанавливаю h равным 10E-14 или ниже, Df начинает получать значения, которые действительно далеки от ожидаемого значения f(10), и относительная ошибка между ожидаемым значением и Df становится огромной.

Почему это?Что здесь происходит?

Answer 1

Оценка f(x) в лучшем случае имеет ошибку округления |f(x)|*mu, где mu - машинная постоянная типа с плавающей запятой.Общая ошибка формулы центральной разности, таким образом, приблизительно равна

2*|f(x)|*mu/(2*h)  +  |f'''(x)|/6 * h^2

. В данном случае экспоненциальная функция равна всем ее производным, так что ошибка пропорциональна

mu/h + h^2/6

*.1009 *, который имеет минимум h = (3*mu)^(1/3), который для двойного формата с mu=1e-16 равен h=1e-5.

Точность увеличивается, если вместо 2*h фактическая разница (x+h)-(x-h) междуоценка баллов используется в знаменателе.Это можно увидеть на следующем графике журнала о расстоянии до точной производной.

Answer 2

В дополнение к ответу @ LutzL Я добавлю некоторую информацию из великой книги Числовые рецепты 3-е издание: Искусство научных вычислений из главы 5.7 о Числовое Производные , особенно о выборе оптимального h значения для данного x:

• Всегда выбирайте h, чтобы h и x отличались точно представимым числом. Следует избегать таких забавных вещей, как 1/3, за исключением случаев, когда x равен чему-то по линии 14.3333333.
• Ошибка округления составляет примерно epsilon * |f(x) * h|, где epsilon - это точность с плавающей запятой, Python представляет числа с плавающей запятой с двойной точностью, поэтому 1e-16. Может отличаться для более сложных функций (где ошибки точности возникают далее), хотя это не ваш случай.
• Выбор оптимального h: Не вдаваясь в подробности, это будет sqrt(epsilon) * x для простого форвардного случая, за исключением случаев, когда ваш x близок к нулю (вы найдете больше информации в книге), что является вашим случаем. Вы можете использовать более высокие значения x в таких случаях, дополнительный ответ уже предоставлен. В случае f(x+h) - f(x-h), как в вашем примере, оно будет составлять epsilon ** 1/3 * x, то есть приблизительно 5e-6 умноженное на x, что может оказаться немного сложным в случае небольших значений, подобных вашему. Совершенно близко (если можно так сказать, имея в виду арифметику с плавающей запятой ...) к практическим результатам, опубликованным @ LutzL хотя.
• Вы можете использовать другие производные формулы, кроме symmetric, который вы используете. Возможно, вы захотите использовать оценку forward или backward (если функция является дорогостоящей для оценки, и вы заранее вычислили f(x). Если ваша функция оценивается дешево, вы можете оценить ее несколько раз, используя более высокий порядок способы уменьшения погрешности точности (см. трафарет из пяти пунктов в википедии , как указано в комментарии к вашему вопросу).

Answer 3

Возможно, вы столкнулись с некоторой численной нестабильностью, например, для x = 10 и h = ~ 1E-13 аргумент для np.exp очень близок к 10, независимо от того, прибавляется или вычитается h, поэтому небольшие ошибки аппроксимации в значении np.exp значительно масштабируется делением с очень маленькими 2 * ч.

Answer 4

В этом руководстве по Python объясняется причина ограниченной точности. Таким образом, десятичные дроби в конечном итоге представлены в двоичном формате, а точность составляет около 17 значащих цифр. Итак, вы правы, что после 10E-14 это становится нечетким.