Численное дифференцирование: аппроксимация производных и конечные разности

Когда необходимо численное дифференцирование?

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

Численное дифференцирование становится необходимым в двух основных случаях. Во-первых, когда функция задана алгоритмически (например, в виде сложной компьютерной подпрограммы или черного ящика), и получить ее производную аналитически невозможно. Во-вторых, когда функция задана в виде дискретного набора точек (таблично) по результатам физического эксперимента, наблюдений или показаний датчиков. В этом случае мы просто не имеем формулы, которую можно было бы продифференцировать. Вычисление скорости объекта по таблице его координат или вычисление ускорения по таблице скоростей — типичные задачи численного дифференцирования.

Формулы конечных разностей

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

• Правая (левая) разность. Вычисляется по формуле: f_prime(x) приближенно равно (f(x+h) - f(x)) / h. Это простейшая формула, имеющая порядок точности O(h). Она несимметрична и вносит определенное смещение в результат вычислений, но часто используется на границах интервалов, где данные с одной из сторон отсутствуют. Аналогично записывается левая разность: (f(x) - f(x-h)) / h.
• Центральная разность. Эта формула использует точки по обе стороны от точки интереса: f_prime(x) приближенно равно (f(x+h) - f(x-h)) / (2*h). За счет симметрии центральная разность обладает гораздо более высокой точностью — второй порядок аппроксимации O(h^2). Это означает, что при уменьшении шага h в 10 раз, ошибка метода уменьшится в 100 раз. Центральные разности являются золотым стандартом в численном моделировании при решении дифференциальных уравнений в частных производных.

Парадокс численного дифференцирования

На первый взгляд кажется, что для получения идеальной точности нужно просто взять шаг h как можно ближе к машинному нулю. Однако здесь возникает знаменитый парадокс численного дифференцирования: операция дифференцирования является некорректно поставленной задачей по Адамару.

При уменьшении шага h погрешность метода (усечения ряда Тейлора) действительно уменьшается. Но в формулах разностей в числителе мы вычитаем два очень близких числа: f(x+h) и f(x). При ограниченной разрядности ЭВМ это приводит к так называемой потере значащих цифр из-за катастрофического взаимного уничтожения. Разность близких чисел обрастает «цифровым шумом», а затем эта ошибка делится на очень малое число h, что приводит к гигантскому взрыву погрешности округления.

Таким образом, график полной погрешности имеет форму буквы «V». Слева погрешность велика из-за ошибок округления, справа — из-за математического усечения. Задача инженера или математика — найти оптимальный шаг h, находящийся на дне этой впадины, где сумма обеих ошибок минимальна. Для 64-битных систем (double precision) оптимальный шаг обычно лежит в районе квадратного корня из машинного эпсилон (примерно 10^-8).