Оптимальный шаг дифференцирования

При дифференцировании функции имеются два основных источника погрешности

• Ошибка метода — уменьшается при уменьшении \(h\)

• Ошибка вычислений — растет при уменьшении \(h\)

Поскольку характер роста ошибок различный, существует некоторое значение \(h^*\), при котором ошибка минимальна. Рассмотим полную ошибку

\[ \varepsilon_\text{total} = \varepsilon_\text{method} + \varepsilon_\text{round} \]

как функцию от \(h\) и найдем минимум.

Схема 1-го порядка

\[ f^{\prime}(x) \approx \frac{f(x+h)-f(x)}{h} \]

Продифферецируем полную ошибку

\[\begin{split} \varepsilon_{\text{total}}(h) = \frac{M_2h}{2} + \frac{2\Delta f}{h}\\ \end{split}\]

по \(h\):

\[\begin{split} 0 = \varepsilon_\text{total}'(h^*) = \frac{M_2}{2} - \frac{2\Delta f}{{h^*}^2}\\h^* = 2\sqrt\frac{\Delta f}{M_2} \end{split}\]

Для функции \(f(x) = \sin x\) оценки производных \(M_2 = M_3 = 1\). Также примем \(\Delta f = 10^{-16}\) - аля работаем с типом double и \(sin x\) можно максимально оценить единицей. Тогда

\[ h^* = 2 \cdot 10^{-8}, \qquad \varepsilon_{\text{total}}^* = 2\cdot 10^{-8} \]

Схема 2-го порядка

\[ f^{\prime}(x) \approx \frac{f(x+h)-f(x-h)}{2 h} \]

Проделав то же самое для формулы дифференцирования второго порядка, получаем

\[\begin{split} \varepsilon_{\text{total}}(h) = \frac{M_3h^2}{6} + \frac{2\Delta f}{2h}\\ \end{split}\]

по \(h\):

\[\begin{split} 0 = \varepsilon_\text{total}'(h^*) = \frac{M_3h^*}{3} - \frac{\Delta f}{{h^*}^2}\\h^* = \sqrt[3]{\frac{3\Delta f}{M_3}} \end{split}\]

При тех же значениях \(M_2, M_3\) и \(\Delta f\) получаем

\[ h^* \approx 6.69 \cdot 10^{-6}, \qquad \varepsilon_{\text{total}}^* \approx 2.24\cdot 10^{-11} \]

Т.е. оптимально брать меньший шаг и при этом будем получать меньшую погрешность!

import numpy as np
import matplotlib.pyplot as plt

def diff1(f, x0, h):
    return (f(x0 + h) - f(x0)) / h

def diff2(f, x0, h):
    return (f(x0 + h) - f(x0 - h)) / (2 * h)

hs = np.logspace(-16, 0, num=50) # h = 1e-16 ... 1
errs1 = []
errs2 = []

for h in hs:
    errs1.append(abs(diff1(np.sin, 1, h) - np.cos(1)))
    errs2.append(abs(diff2(np.sin, 1, h) - np.cos(1)))

M2 = M3 = 1
plt.figure(figsize=(10, 7))
plt.loglog(hs, errs1, 'k.', ms=16)
plt.loglog(hs, errs2, 'r.', ms=16)
plt.loglog(hs, M2 * hs / 2, 'k-', label='First order', lw=3)
plt.loglog(hs, M3 * hs**2 / 6, 'r-', label='Second order', lw=3)
plt.ylim(1e-14, 1)
plt.xlabel('$h$')
plt.ylabel('$\\varepsilon$')
plt.legend(loc='upper center')
plt.show()