Численное дифференцирование: конечные разности и выбор шага

Датчик пишет скорость насоса раз в секунду, производной в этой записи нет — есть столбик чисел. Или хуже: функция существует только в виде программы, и аналитическая формула недостижима в принципе. Производную придётся выжимать из значений. Парадокс в том, что определение через предел — $h\to 0$ — здесь не помощник, а вредитель: чем меньше шаг, тем ближе мы к машине, которая считает с конечной точностью. Найдём золотую середину.

Определение через предел и его численная ловушка#

По определению $f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}$. Подставим $h=10^{-15}$ и попробуем найти производную $\sin x$ в точке $x=1$. Истинный ответ $\cos 1=0.5403023$. А машина вернёт $0.5551$ — и это ещё не худший случай: при $h=10^{-16}$ ответ ровно ноль.

Разбор полёта. В double числа около единицы хранятся с точностью $\epsilon \approx 2.2\cdot 10^{-16}$, значит, $\sin (1+10^{-16})$ — соседнее машинное число, синус которого неотличим от $\sin 1$. Числитель обнуляется. При $h=10^{-15}$ числитель не ноль, но состоит почти целиком из шума округления: два близких значения функции вычитаются, старшие цифры гасят друг друга — то самое катастрофическое сокращение из урока про источники погрешностей. Вывод: определение в лоб неприменимо, нужен компромиссный шаг.

Разности вперёд, назад и центральная#

Ряд Тейлора раскладывает всё по полочкам: $f(x+h)=f(x)+h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)+\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+\dots$. Каждый шаблон разностей — это способ так скомбинировать значения функции, чтобы лишние члены сократились. Сдвинем ряд на $-h$ и вычтем его из исходного: члены с $f$ и $f^{\prime \prime }$ взаимно уничтожатся, и останется производная с добавкой порядка $h^2$.

Разность назад $\frac{f(x)-f(x-h)}{h}$ симметрична вперёд и тоже даёт $O(h)$. Почему центральная вдвое точнее? Симметрия: сдвиги $+h$ и $-h$ вносят поправки с противоположными знаками, и все чётные степени шага выпадают. Первый неисправленный член — $\frac{h^2}{6}{f}^{\prime \prime \prime }$.

1. 1
   Пишем оба разложения: $f(x+h)=f+h{f}^{\prime }+\frac{h^2}{2}{f}^{\prime \prime }+\frac{h^3}{6}{f}^{\prime \prime \prime }+\dots$ и $f(x-h)=f-h{f}^{\prime }+\frac{h^2}{2}{f}^{\prime \prime }-\frac{h^3}{6}{f}^{\prime \prime \prime }+\dots$
2. 2
   Вычитаем второе из первого: члены с $f$ и $f^{\prime \prime }$ сокращаются, а нечётные степени удваиваются: $f(x+h)-f(x-h)=2h{f}^{\prime }+\frac{h^3}{3}{f}^{\prime \prime \prime }+\dots$
3. 3
   Делим на $2h$: $f^{\prime }+\frac{h^2}{6}{f}^{\prime \prime \prime }+\dots$ — чётные степени шага исчезли, точность $O(h^2)$
4. 4
   Контроль на числах: для $\sin x$ в $x=1$ при $h=0.1$ формула даёт $0.539402$ против $\cos 1=0.540302$ — ошибка $9\cdot 10^{-4}$, ровно $h^2\cdot |\cos 1|/6=0.0009$

Пять точек вместо трёх: шаблон повышенного порядка

Если функция доступна в обе стороны и на два шага, точность можно поднять до $O(h^4)$ без уменьшения шага — комбинируя пять значений с весами, гасящими и $h^2$, и $h^3$ члены:

Проверка на $\sin x$ в $x=1$ при $h=0.1$: пятиточечная схема даёт $0.540300507$ — ошибка $1.8\cdot 10^{-6}$ против $9\cdot 10^{-4}$ у трёхточечной на том же шаге: в пятьсот раз точнее за те же два лишних вызова функции. Но пятиточечный шаблон унаследовал и главный порок: деление на $h$ усиливает шум так же, а весов стало больше. Уменьшать шаг до бесконечности всё равно нельзя — теперь за дело берётся следующая глава.

Вторые производные: центральная схема и её цена#

Сложим ряды для $f(x+h)$ и $f(x-h)$ и вычтем удвоенное $f(x)$: первые степени сократятся, вторые — удвоятся. Деление на $h^2$ даёт шаблон для второй производной:

Проверим на $\sin x$ в $x=1$: истинная вторая производная $-\sin 1=-0.841470985$. При $h=0.1$ схема даёт $-0.840770$ (ошибка $7\cdot 10^{-4}$), при $h=0.01$ — $-0.841463972$ (ошибка $7\cdot 10^{-6}$): шаг в десять раз меньше, ошибка в сто раз меньше, порядок $O(h^2)$ подтверждён. Но у этой схемы порок старшего брата: деление на $h^2$ усиливает шум округления вдвое быстрее, и оптимум шага сдвигается к $h\approx 10^{-4}$ — дальше ошибка только растёт. При $h=10^{-8}$ ответ $-1.11$ вместо $-0.8415$: мусор победил.

Дилемма шага: методическая ошибка против округления#

Вот сердце всей темы. Методическая ошибка центральной разности равна $\frac{h^2}{6}|f^{\prime \prime \prime }|$ и падает с уменьшением шага. Ошибка округления ведёт себя противоположно: вычитание близких чисел оставляет шум порядка $\epsilon |f|$, а деление на малое $h$ его раздувает — примерно $\frac{\epsilon |f|}{h}$. Сумма двух конкурентов имеет минимум в точке, где они сравниваются: $h\sim \sqrt[3]{\epsilon }\approx 6\cdot 10^{-6}$ для double.

Таблица снята с настоящей арифметики double, и кривая ошибки V-образна: слева методика, справа шум, дно — у $h\approx 10^{-5}$. Односторонняя разность живёт по той же схеме, но её оптимум сдвинут к $h\approx \sqrt{\epsilon }\approx 10^{-8}$, а потолок точности — порядка $10^{-9}$ против $10^{-11}$ у центральной. Дороже ли платить за точность? Центральной схеме нужны значения с двух сторон — на краях таблицы их просто нет, и приходится довольствоваться односторонними шаблонами.

Производная по таблице данных#

Таблица с постоянным шагом $h$ развязывает руки: шаблоны переносятся на узлы буквально. Для узла $i$ центральная разность — $\frac{y_{i+1}-y_{i-1}}{2h}$; на левом крае таблицы нет $y_{i-1}$, и работает разность вперёд $\frac{y_{i+1}-y_i}{h}$, на правом — назад. Неравномерный шаг усложняет коэффициенты, но не меняет идею: комбинируем соседние значения так, чтобы члены Тейлора сокращались.

Числа для ощутимости. Пусть $y_i=\sin (x_i)$ на сетке с шагом $0.2$, и нужна производная в узле $x=0.8$. Истина $\cos 0.8=0.6967$. Разность вперёд $\frac{\sin 1-\sin 0.8}{0.2}=0.6206$ — ошибка $7.6\cdot 10^{-2}$. Центральная $\frac{\sin 1-\sin 0.6}{0.4}=0.6921$ — ошибка $4.6\cdot 10^{-3}$, в шестнадцать раз лучше на тех же данных. Один и тот же столбик чисел, одна и та же цена — разная точность.

• Внутренние точки — центральная разность: она даёт $O(h^2)$ без всяких дополнительных измерений
• Края таблицы — односторонние шаблоны, а точнее их версии второго порядка: $\frac{-3f(x)+4f(x+h)-f(x+2h)}{2h}$ тоже даёт $O(h^2)$
• Шумные данные — сначала сглаживание или МНК, потом разности: дифференцирование усиливает шум как $\delta /h$
• Нужна вторая производная — помните про $h\approx 10^{-4}$ как оптимум и не верьте цифрам при $h\le 10^{-6}$
• Проверка всегда: сравните численный ответ с аналитическим эталоном на тестовой функции, где производная известна

Отсюда производственный рецепт: если таблица зашумлена — строим по ней гладкую модель и дифференцируем модель. Кандидат на роль модели — интерполяционный сплайн по точкам или регрессия из урока про метод наименьших квадратов. Обратная операция — интегрирование — шум, наоборот, глушит; этому посвящён урок численное интегрирование. Определение производной и её геометрию напомнит урок производная: определение и касательная, из которого всё и выросло; разложения, на которых держатся выводы порядков, разобраны в формуле Тейлора, а для вторых производных высших порядков пригодится материал о производных высших порядков. Перед сдачей работы сверьте таблицу производных в шпаргалке — численный ответ обязан ложиться на аналитический эталон.