Numerical differentiation. Calculation of the first derivative.¶
Definition¶
By definition, the first derivative of a smooth function \(f(x)\) at a point x is calculated as
When calculating the first derivative of the function \(f(x)\) on a computer, we replace the infinitesimal \(h \rightarrow \infty\) with a small but finite value \(h\):
where \(\mathrm{O}(h)\) is the derivative calculation error, which naturally depends on \(h\). Previous formula is called a difference scheme for calculating the first derivative (more precisely, a right difference scheme or just a right difference). Similarly, maybe the left-hand difference scheme is written.
How to determine \(\mathrm{O}(h)\)? Expand the function \(f(x)\) in a Taylor series at the point \(x + h\):
whence it follows that in the first order of the expansion
By choosing a very small \(h\), the round-off errors in computing on a computer can be comparable to or greater than \(h\). Therefore, we are interested in an algorithm that gives lower error value for the same value of \(h\).
Such an improved algorithm can be easily obtained by expanding the function \(f(x)\) into a Taylor series at the points \(x + h\) and \(x - h\), then subtracting one result from the other, which gives
where the error in calculating the first derivative
This is the central difference scheme (central difference).
In principle, it is possible to follow the path of improving the accuracy of the method for calculating the first derivative and further. For example, considering the expansion of the function \(f(x)\) in a Taylor series at the points \(x + h\), \(x + 2h\), \(x - h\), and \(x - 2h\), one can obtain a four-point scheme etc.
Usage¶
Imagine that we want to find the derivative of the following function:
Then the code will look like this:
// 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;
}