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

Question

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

Я реализую метод Ньютона-Рафсона , и он требует получения производной функции.

Answer 1

Я согласен с @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.

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

Answer 2

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

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

Answer 3

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

- относительная погрешность округления, около2 ^ {- 16} для двойного и 2 ^ {- 7} для числа с плавающей запятой.

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

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

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

Обратите внимание, что ошибка второго приближения стремится к 0 быстрее, чем первая.Но если f '' '(x) очень запаздывает, то первый вариант более предпочтителен:

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

Answer 4

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

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

Answer 5

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

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

Answer 6

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

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

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

Answer 7

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

Answer 8

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

Пол