Численное дифференцирование

[andreso]

Вопрос абстрактный, но по другому его не задашь.
Имеется алгоритмически заданная функция. Функция гарантированно гладкая.
Нужно найти ее производную численно. При использовании метода центральных разностей столкнулся с неопределенностью необходимого шага дифференцирования. Т.е. фактически при различных шагах получаются различные результаты. Выбрать шаг нужно автоматически, поскольку при расчетах параметры функции могут меняться. Как это лучше сделать?

[ewert]

andreso в сообщении #748816 писал(а):

Выбрать шаг нужно автоматически, поскольку при расчетах параметры функции могут меняться. Как это лучше сделать?

Полностью автоматически -- никак. Но практически приемлемые результаты получатся, если взять для начала шаг, допустим, $10^{-3}$, потом пересчитывать значение производной, последовательно уменьшая шаг вдвое (например) -- до тех пор, пока разности между соседними итерациями не перестанут уменьшаться по величине.

Значение $10^{-3}$ взято, естественно, с потолка; просто оно характерно. В типичных случаях оптимум достигается при $h\sim10^{-5}...10^{-6}$ (т.е. понадобится порядка 7-10 итераций), а соответствующая погрешность оказывается порядка $h\sim10^{-10}...10^{-11}$. Большей точности (в стандартных режимах вычислений) с помощью этой формулы получить невозможно.

[VPro]

Лучше использовать известную трехточечную формулу:
$f'_i = \frac{-3 f_{i-1} + 4 f_i - f_{i+1}}{2 h}$
Ошибка вычисления равна $ \frac{f'''(\xi)}{3} h^2$, где $\xi$ расположена между $i-1$-тым и $i+1$ узлами. Ограничив значение ошибки и оценив сверху 3-ю производную, Вы получите автоматический критерий выбора шага.

[ewert]

VPro в сообщении #748831 писал(а):

Лучше использовать известную трехточечную формулу:
$f'_i = \frac{-3 f_i + 4 f_{i+1} - f_{i+2}}{2 h}$

Ничего хорошего в этой формуле нет, кроме того, что она при том же порядке менее точна чем та, симметричная (которая, кстати, формально тоже трёхточечная). Если же хочется повысить точность, то надо использовать, скажем, симметричную четырёхточечную (формально пятиточечную) формулу $f'_0 = \frac{f_{-2}-8 f_{-1} + 8 f_{1} - f_{2}}{12 h^2}+O(h^4)$.

VPro в сообщении #748831 писал(а):

Ограничив значение ошибки и оценив сверху 3-ю производную, Вы получите автоматический критерий выбора шага.

Ничего Вы не получите, т.к. отсутствует оценка погрешностей округления, а именно они здесь критичны. Не говоря уж о том, что оценивать третью производную в условиях, когда нет ещё даже первой -- это любопытно.

[VPro]

1) С первой частью согласен. Можно с равным успехом использовать и ту и другую формулы. Я привык к трехточечной. И, бесспорно, пятиточечная формула еще лучше.

2) По прежнему полагаю, что автоматический выбор шага вполне реален. Пусть, например, $f'''<10^3$. Тогда, чтобы ошибка вычисления производной не превосходила $10^{-2}$ (для производной вполне достаточно), нужно взять шаг порядка $10^{-3}$. С другой стороны, погрешность округления, если использовать 8-байтовой представление чисел, дает 15 верных знаков, а 10-байтовое дает 19 верных знаков. Не вижу критичности этой погрешности. Впрочем, конечно, все это зависит от конкретной функции и требований к точности.

[ewert]

VPro в сообщении #748847 писал(а):

Я привык к трехточечной.

Формула $f'_i = \frac{f_{i+1}- f_{i-1}}{2 h}+O(h^2)$ -- тоже трёхточечная. Однако она проще (что малосущественно) и при этом несколько точнее (что уже может оказаться полезным, хоть и не сильно)

VPro в сообщении #748847 писал(а):

10-байтовое дает 19 верных знаков. Не вижу критичности этой погрешности.

Вообще-то 17-18; и ещё собственно вычисление функции повышает погрешность.

Критично или нет -- зависит от задачи. Это совсем не существенно в разных сеточных схемах, где шаг заведомо не слишком мал просто по смыслу алгоритма. Если же нужны по возможности точные значения производной, то тут дело другое. Во всяком случае, полезно иметь в виду, что даже в десятибайтовом случае невозможно получить трёхточечными формулами более чем 12 значащих цифр.

[VPro]

ewert в сообщении #748852 писал(а):

...
Если же нужны по возможности точные значения производной, то тут дело другое. Во всяком случае, полезно иметь в виду, что даже в десятибайтовом случае невозможно получить трёхточечными формулами более чем 12 значащих цифр.

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

-- Ср июл 24, 2013 14:26:14 --

Да, и Вы меня убедили в преимуществе формулы центральной разности. Буду использовать ее.

[worm2]

ewert в сообщении #748832 писал(а):

Ничего Вы не получите, т.к. отсутствует оценка погрешностей округления...

Ну и чтобы окончательно сбить с толку топикстартера добавлю, что поскольку неизвестно, как функция вычисляется, то неизвестна и та часть её погрешности, которая обусловлена округлениями, накапливающимися в ходе её вычисления (про другие части погрешности я даже и не заикаюсь).

Поэтому можно ожидать (ну, если не возникнет ещё каких-нибудь проблем ), что при уменьшении шага численно посчитанная производная будет сначала "как бы сходиться" к какому-то значению, а потом (заранее нельзя сказать когда), вдруг начнёт "плясать", сначала слабо, а потом всё сильнее, затем может впасть в ноль, либо уверенно увеличиваться по модулю (возможен также сразу резкий спад до нуля, без плясок), пока шаг не станет меньше минимально возможного при выбранном формате хранения чисел (этот минимум, к счастью, легко вычисляется), после чего случится деление на ноль. Если вспомнить, что разум — это вялотекущая шизофрения, аналогия получается полной

Наша задача — поймать момент, пока производная не начала "плясать" сильнее, чем "сходиться", это и будет самым точным значением. Практическим критерием иногда можно считать локальный экстремум (смену убывания на возрастание или наоборот). Но тут есть нюансы. Например, если производная близка к нулю, нужно как-то уметь отличать внезапное падение в ноль от плавного уменьшения до нуля. В первом случае мы в качестве производной берём последнее ненулевое значение, во втором — 0.

[Portnov]

Вообще-то, численное дифференцирование — некорректно поставленная (по Адамару) задача. Практические результаты этого указал worm2: маленькие погрешности в вычислении исходной функции могут привести к сколь угодно большим погрешностям при вычислении производной.

[ewert]

worm2 в сообщении #749053 писал(а):

Наша задача — поймать момент, пока производная не начала "плясать" сильнее, чем "сходиться", это и будет самым точным значением. Практическим критерием иногда можно считать локальный экстремум (смену убывания на возрастание или наоборот). Но тут есть нюансы.

Тут я бы сказал два ньюанеца.

Один состоит в том, что значение разностной производной может на какое-то время стабилизироваться на значении, равном истинной производной (в пределах погрешности округления). Случайно так вполне может произойти, и происходит достаточно часто. Это сильно препятствует оценке достигнутой точности, но никак не нарушает сам алгоритм поиска оптимума.

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

[VinSer100]

http://fourier.eng.hmc.edu/e176/lectures/ch6/node9.html

http://ceres-solver.org/numerical_derivatives.html

[ilghiz]

andreso в сообщении #748816 писал(а):

Вопрос абстрактный, но по другому его не задашь.
Имеется алгоритмически заданная функция. Функция гарантированно гладкая.
Нужно найти ее производную численно. При использовании метода центральных разностей столкнулся с неопределенностью необходимого шага дифференцирования. Т.е. фактически при различных шагах получаются различные результаты. Выбрать шаг нужно автоматически, поскольку при расчетах параметры функции могут меняться. Как это лучше сделать?

Если число арифметических операций при вычислении этой функции не приближается к числу слов с плавающей точкой, которые можно дополнительно выделить в оперативной памяти, я бы в любом случае, попользовал метод Баура-Штрассена. Его арифметическая сложность только в разы больше, чем сложность вычисления самой функции, можно написать короткий С++ класс и перегрузить все операторы.

[novichok2018]

Попробуйте метод двойного пересчёта Рунге. Там и погрешность примерно оценивается.