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

Я согласен с @erikkallen в том, что (f(x + h) - f(x - h)) / 2 * h - это обычный подход для численного приближения производных. Однако получить правильный размер шага h немного сложно.

Ошибка аппроксимации в (f(x + h) - f(x - h)) / 2 * h уменьшается по мере того, как h становится меньше, что говорит о том, что вы должны брать h как можно меньше. Но когда h становится меньше, ошибка вычитания с плавающей запятой увеличивается, поскольку числитель требует вычитания почти равных чисел. Если h слишком мало , вы можете потерять большую точность при вычитании. Поэтому на практике вам нужно выбрать не слишком маленькое значение h, которое сводит к минимуму сочетание ошибки аппроксимации и числового ошибка.

Как правило, вы можете попробовать h = SQRT(DBL_EPSILON), где DBL_EPSILON — наименьшее число с двойной точностью e, такое что 1 + e != 1 с машинной точностью. DBL_EPSILON примерно соответствует 10^-15, поэтому вы можете использовать h = 10^-7 или 10^-8.

Дополнительные сведения см. в этих примечаниях по выбору размера шага для дифференциальных уравнений.

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

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

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

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

Мы можем рассчитать общую ошибку:

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

2) Второй случай:

Обратите внимание, что ошибка второго приближения стремится к 0 быстрее, чем ошибка первого. Но если 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''(x) или f'''(x).

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

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

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

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

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

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

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

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

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

Павел