Численное дифференцирование. Расчет первой производной.

Определение

По определению первая производная гладкой функции \(f\) вточке x вычисляется как

\begin{equation} f^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h}. \end{equation}

При вычислении первой производной функции \(f(x)\) на компьютере мы заменяем бесконечно малую \(h \rightarrow \infty\):

\begin{equation} f^{\prime}(x)=\frac{f(x+h)-f(x)}{h}+\mathrm{O}(h), \end{equation}

где \(\mathrm{O}(h)\) - это ошибка вычисления производной, которая, естественно, зависит от \(h\). Предыдущая формула называется разностной схемой для вычисления первой производной (точнее, правильной разностной схемой или просто правильной разностной). Точно так же, может быть, написана левая разностная схема.

Как определить \(\mathrm{O}(h)\)? Разложим функцию \(f(x)\) в ряд Тейлора в точке \(x + h\):

\begin{equation} 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)+\ldots, \end{equation}

откуда следует, что в первом порядке разложения

\begin{equation} \mathrm{O}(h)=-\frac{h}{2} f^{\prime \prime}(x)+\ldots. \end{equation}

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

Такой улучшенный алгоритм может быть легко получен путем разложения функции \(f(x)\) в ряд Тейлора в точках \(x + h\) и \(x - h\), а затем вычитания одного результата из другого, что дает

\begin{equation} f^{\prime}(x)=\frac{f(x+h)-f(x-h)}{2 h}+\mathrm{O}\left(h^{2}\right), \end{equation}

где погрешность вычисления первой производной

\begin{align*} \mathrm{O}\left(h^{2}\right)=-\frac{h^{2}}{6} f^{\prime \prime \prime}(x)+\ldots. \end{align*}

Это центральная разностная схема (центральная разностная схема).

В принципе, можно пойти по пути повышения точности метода вычисления первой производной и далее. Например, рассматривая разложение функции f (x) в ряд Тейлора в точках \(x + h\), \(x + 2h\), \(x - h\) и \(x - 2h\), можно получить четырехточечную схему и т. д.

Использование

Представьте, что мы хотим найти производную следующей функции:

\begin{equation} f(x) = \sin{(x)} \end{equation}

Тогда код будет выглядеть так:

// example_first_order_derivative_h.cpp

#include <iostream>
#include "../src/numerary.hpp" // Numerary library

using namespace std;
using namespace numerary;

/* Functiion to derive */
double f(double x) {
    return sin(x);
}

/* The main function */
int main() {

    const short int order = 1;
    double x, dy_dx;

    // Point where we want get value of derivative function
    x = M_PI;

    dy_dx = Numerary::differentiate(f, order, x);

    cout << "dy/dx (" << x << ") = " << dy_dx << endl;

    return 0;
}