Конечные разности: производная по точкам

Урок про численное дифференцирование: как приблизить производную разностью значений и почему здесь нельзя «просто взять h поменьше».

Конечная разность приближает производную $f'(x)$ отношением приращения функции к приращению аргумента, например $\frac{f(x+h) - f(x)}{h}$.

Идея: из определения производной

Производная по определению — предел $\frac{f(x+h) - f(x)}{h}$ при $h \to 0$. Численно мы не можем взять предел, но можем взять конкретное малое h — получится правая (односторонняя) разность. Она проста, но не очень точна. Гораздо лучше центральная разность $\frac{f(x+h) - f(x-h)}{2h}$: симметрия гасит главную часть ошибки, и точность подскакивает с первого порядка до второго.

import math

f = lambda x: math.sin(x)
x0 = 1.0
истина = math.cos(x0)        # производная sin = cos

print(f"{'h':>8} | {'правая':>12} | {'центральная':>14}")
print("-" * 40)
for h in [0.1, 0.01, 0.001, 1e-5]:
    правая = (f(x0 + h) - f(x0)) / h
    центр  = (f(x0 + h) - f(x0 - h)) / (2 * h)
    print(f"{h:8g} | {abs(правая - истина):12.2e} | {abs(центр - истина):14.2e}")

Вывод:

       h |       правая |    центральная
----------------------------------------
     0.1 |     4.29e-02 |       9.00e-04
    0.01 |     4.22e-03 |       9.00e-06
   0.001 |     4.21e-04 |       9.01e-08
   1e-05 |     4.21e-06 |       1.11e-11

Сравните столбцы. У правой разности при делении h на 10 ошибка падает в 10 раз — это первый порядок ($\text{ошибка} \propto h$). У центральной при делении h на 10 ошибка падает в 100 раз — второй порядок ($\text{ошибка} \propto h^2$). При $h = 0.001$ центральная разность уже на 4 порядка точнее правой при той же стоимости.

Порядок точности и его источник

Откуда берутся эти порядки? Из ряда Тейлора. Разложив $f(x+h)$, видим, что в правой разности остаётся член $\sim \frac{h}{2} \cdot f''$ — отсюда первый порядок. В центральной разности симметричные члены с $f''$ сокращаются, и первым ненулевым остатком оказывается $\sim \frac{h^2}{6} \cdot f'''$ — отсюда второй порядок. Этот же приём (комбинировать значения в нескольких точках, чтобы сократить младшие члены Тейлора) даёт формулы ещё выше: пятиточечная схема имеет четвёртый порядок.

Почему дифференцирование «шумит»

Казалось бы, бери h совсем крошечным — и точность бесконечна. Но нет: численное дифференцирование — неустойчивая операция. При малом h мы вычитаем два почти равных значения $f(x+h)$ и $f(x-h)$ (потеря значимости!) и делим на крошечное h — ошибка округления взрывается как $\sim \frac{\varepsilon}{h}$. Складываясь с погрешностью метода $\sim h^2$, она даёт классическую «долину»: есть оптимальный h, после которого точность не растёт, а падает.

import math
истина = math.e            # производная e^x в точке 1 = e

for k in [2, 4, 6, 8, 10, 12]:
    h = 10.0 ** (-k)
    производная = (math.exp(1.0 + h) - math.exp(1.0 - h)) / (2 * h)
    print(f"h=1e-{k:<2}: ошибка = {abs(производная - истина):.2e}")

Вывод:

h=1e-2 : ошибка = 4.53e-05
h=1e-4 : ошибка = 4.53e-09
h=1e-6 : ошибка = 1.63e-10
h=1e-8 : ошибка = 6.60e-09
h=1e-10: ошибка = 6.73e-07
h=1e-12: ошибка = 2.10e-04

Ошибка падает до $h \approx 10^{-6}$ (минимум $\sim 1.6 \times 10^{-10}$), а затем растёт: при $h = 10^{-12}$ она хуже, чем при $h = 10^{-2}$! Это прямое следствие потери значимости. Оптимальный h для центральной разности — порядка $\varepsilon^{1/3} \approx 10^{-5} \ldots 10^{-6}$, а не «как можно меньше».

Как работает под капотом

Полная ошибка центральной разности $\approx C_1 \cdot h^2 + C_2 \cdot \frac{\varepsilon}{h}$. Первый член (метода) убывает с h, второй (округления) растёт. Минимум суммы достигается, когда они сравнимы — отсюда оценка оптимального $h \sim \varepsilon^{1/3}$ для второго порядка. Практический вывод: для производных не гонитесь за микроскопическим шагом. А если функция задана зашумлёнными данными (не формулой), численное дифференцирование вообще опасно — оно усиливает шум; там сначала сглаживают (например, аппроксимируют МНК и дифференцируют гладкую модель). Современная альтернатива для гладких функций, заданных кодом, — автоматическое дифференцирование, дающее производную без шага h и без ошибки усечения.

Частые ошибки

Брать h как можно меньше. Ниже оптимума округление губит точность; для центральной разности оптимум $\sim 10^{-5} \ldots 10^{-6}$.
Дифференцировать зашумлённые данные напрямую. Разности усиливают шум; сначала сглаживайте.
Использовать одностороннюю разность без нужды. Центральная вдвое-вчетверо точнее при той же цене (если доступны обе стороны).

Итоги

Конечные разности приближают производную отношением приращений; центральная (2-й порядок) точнее правой (1-й).
Порядок выводится из ряда Тейлора сокращением младших членов.
Численное дифференцирование неустойчиво: при малом h округление ($\sim \frac{\varepsilon}{h}$) перевешивает метод ($\sim h^2$).
Существует оптимальный h ($\sim \varepsilon^{1/3}$); «меньше» $\ne$ «точнее».