Реализация производной в C/С++
Как производная от f(x)
обычно рассчитывается программно для обеспечения максимальной точности?
Я реализую метод Newton-Raphson и требует принятия производной функции.
Ответы
Ответ 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.
Подробнее см. в этих примечаниях о выборе размера шага для дифференциальных уравнений.
Ответ 2
Newton_Raphson предполагает, что вы можете иметь две функции f (x) и ее производную f '(x). Если у вас нет функции, доступной как функция, и вам нужно оценить производную от исходной функции, тогда вы должны использовать другой алгоритм поиска корней.
Корневой поиск в Википедии дает несколько предложений, как и любой текст численного анализа.
Ответ 3
![alt text]()
![alt text]()
1) Первый случай:
![alt text]()
- относительная ошибка округления, около 2 ^ {- 16} для double и 2 ^ {- 7} для float.
Мы можем вычислить общую ошибку:
![alt text]()
Предположим, что вы используете операцию двойного плавания. Таким образом, оптимальное значение h равно 2sqrt (DBL_EPSILON/f '' (x)). Вы не знаете f '' (x). Но вы должны оценить это значение. Например, если f '' (x) составляет около 1, то оптимальное значение h равно 2 ^ {- 7}, но если f '' (x) составляет около 10 ^ 6, то оптимальное значение h равно 2 ^ {- 10}!
2) Второй случай:
![alt text]()
Обратите внимание, что ошибка второго приближения стремится к 0 быстрее первого.
Но если f '' '(x) очень лаг, то предпочтительнее первый вариант:
![alt text]()
Заметим, что в первом случае 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).
Ответ 4
fprime(x) = (f(x+dx) - f(x-dx)) / (2*dx)
для некоторого малого dx.
Ответ 5
Что вы знаете о f (x)? Если у вас есть только черный ящик, единственное, что вы можете сделать, это численное приближение производной. Но точность обычно не так хороша.
Вы можете сделать гораздо лучше, если вы можете коснуться кода, который вычисляет f. Попробуйте "автоматическое дифференцирование" . Там есть несколько хороших библиотек для этого. С небольшим количеством магии библиотеки вы можете легко преобразовать свою функцию в то, что автоматически вычисляет производную. Для простого примера С++ см. Исходный код в этом немецком обсуждении.
Ответ 6
Вы определенно хотите учесть предложение Джона Кука для выбора h, но вы, как правило, не хотите использовать разницу по центру, чтобы приблизить производную. Основная причина заключается в том, что она требует дополнительной оценки функции, если вы используете разницу в прямом направлении, то есть
f'(x) = (f(x+h) - f(x))/h
Тогда вы получите значение f (x) бесплатно, потому что вам нужно вычислить его уже для метода Ньютона. Это не так уж и важно, когда у вас есть скалярное уравнение, но если x - вектор, то f '(x) является матрицей (якобиан), и вам нужно будет сделать n дополнительных оценок функций для ее приближения используя подход с разнесением по центру.
Ответ 7
В дополнение к ответу Джона Д. Кукса выше, важно не только учитывать точность с плавающей точкой, но и надежность функции f (x). Например, в финансах обычно случается, что f (x) на самом деле является методом Монте-Карло, а значение f (x) имеет некоторый шум. Используя очень небольшой размер шага, в этих случаях сильно ухудшается точность производной.
Ответ 8
Обычно шум сигнала влияет на качество производных больше, чем что-либо еще. Если у вас есть шум в вашем f (x), Savtizky-Golay - отличный алгоритм сглаживания, который часто используется для вычисления хороших производных. В двух словах SG помещает многочлен локально в ваши данные, тогда этот многочлен можно использовать для вычисления производной.
Пол