Внедрение производной в C/C++

Question

Как производная от f(x) обычно рассчитывается программно для обеспечения максимальной точности?

Я реализую метод Newton-Raphson и требует принятия производной функции.

Answer 1

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

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

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

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

Answer 2

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

для некоторых небольших dx.

Answer 3

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

Wikipedia root finding дает несколько предложений, как и любой текст численного анализа.

Answer 4

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

Вы можете сделать много лучше, если вы можете коснуться кода, который вычисляет f. Попробуйте "automatic differentiation". Там есть несколько хороших библиотек для этого. С небольшим количеством магии библиотеки вы можете легко преобразовать свою функцию в то, что автоматически вычисляет производную. Простой пример C++ см. В немецком обсуждении source code in this.

Answer 5

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

Answer 6

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

- относительная погрешность округления, около 2^{- 16} для двойной и 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 быстрее, чем первый. Но если е «» '(х) очень Lagre, то первый вариант более предпочтителен:

Обратите внимание, что в первом случае ч пропорциональна е, но во втором случае ч пропорциональна е^{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).

Answer 7

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

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

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

Answer 8

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

Paul