Реализация производной в C/C++

Newton_Raphson предполагает, что у вас может быть две функции f (x) и ее производная f '(x). Если у вас нет производной, доступной как функция, и вам необходимо оценить производную от исходной функции, тогда вам следует использовать другой алгоритм поиска корня.

Поиск корня в Википедии дает несколько предложений, как и любой текст численного анализа.

1) Первый случай:

- относительная ошибка округления, около 2 ^ {- 16} для double и 2 ^ {- 7} для float.

Мы можем вычислить общая ошибка:

Предположим, вы используете двойную операцию с плавающей запятой. Таким образом, оптимальное значение h составляет 2sqrt (DBL_EPSILON / f '' (x) ). Вы не знаете f '' (x) . Но вы должны оценить эту стоимость. Например, если f '' (x) примерно равно 1, то оптимальное значение h равно 2 ^ {- 7}, но если f '' (x) составляет примерно 10 ^ 6, тогда оптимальное значение h равно 2 ^ {- 10}! Но если f '' '(x) очень запаздывающее, то первый вариант более предпочтителен:

Обратите внимание, что в первом случае h пропорционально e, а во втором случае h пропорционально e ^ {1/3}. Для двойных операций с плавающей запятой e ^ {1/3} равно 2 ^ {- 5} или 2 ^ {- 6}. (Я полагаю, что f '' '(x) около 1).

Какой способ лучше? Неизвестно, если вы не знаете f '' (x) и f '' '(x) или не можете оценить эти значения. Считается, что второй вариант предпочтительнее. Но если вы знаете, что f '' '(x) очень большое значение, используйте сначала.

Какое оптимальное значение h? Предположим, что f' '(x) и f' '' (x) равны около 1. Также предположим, что мы используем двойные операции с плавающей запятой. Тогда в первом случае h примерно равно 2 ^ {- 8}, в первом случае h примерно равно 2 ^ {- 5}. Исправьте эти значения, если знаете f ''

fprime(x) = (f(x+dx) - f(x-dx)) / (2*dx)

для небольшого dx.

Что вы знаете о f (x)? Если у вас есть только f в виде черного ящика, единственное, что вы можете сделать, - это численно аппроксимировать производную. Но точность обычно не так хороша.

Вы можете сделать намного лучше, если прикоснетесь к коду, который вычисляет f. Попробуйте «автоматическое дифференцирование» . Для этого есть несколько хороших библиотек. Немного библиотечной магии вы можете легко преобразовать свою функцию во что-то, что вычисляет производную автоматически. Простой пример C ++ можно найти в исходном коде в этом немецком обсуждении.

В дополнение к ответу Джона Д. Кука выше важно не только учитывают точность с плавающей запятой, а также надежность функции f (x). Например, в финансах часто бывает, что f (x) фактически является моделированием Монте-Карло, а значение f (x) имеет некоторый шум. Использование очень маленького шага в этих случаях может сильно снизить точность производной.

Вы определенно хотите принять во внимание предложение Джона Кука о выборе h, но обычно вы не хотите использовать центрированную разницу для аппроксимации производной. Основная причина в том, что это требует дополнительной оценки функции, если вы используете прямую разницу, то есть

f'(x) = (f(x+h) - f(x))/h

Тогда вы получите значение f (x) бесплатно, потому что вам нужно вычислить его уже для метода Ньютона. Это не такая уж большая проблема, когда у вас есть скалярное уравнение, но если x - вектор, тогда f '(x) - это матрица (якобиан), и вам нужно будет выполнить n дополнительных вычислений функций, чтобы аппроксимировать ее. используя подход центрированной разности.

Обычно шум сигнала влияет на качество производной больше, чем что-либо еще. Если у вас есть шум в вашем f (x), Савтицкий-Голей - отличный алгоритм сглаживания, который часто используется для вычисления хороших производных. Вкратце, SG подбирает полином локально к вашим данным, затем этот полином можно использовать для вычисления производной.

Пол