Functionals and numerical differentiation

This web page explains one of the mechanisms that you can use in C++ to pass a function as an argument as another function. It also gives as an example a function that attempts to numerically differentiate a function. Some of the shortcomings of this method are discussed. In general differentiation is a very difficult process to carry out numerically.

1. Functionals in C++

This section introduces you to one of the ways you can pass funtion objects to functions in C++. This requires the package <functional> and the type std::function<...> defined there. The following example is based on one from the wikipedia article.

The details require some care. The <functional> package loads the std::function<...> type that allows you to declare a variable as being a function. This was introduced with C++ 2011, and for older compilers it is necessary to change a setting to use 2011 support.

For example, in C++ 2011, a function $f$ that takes a double input and returns a double declared by double f(double x) { ... } and would fit the type std::function<double(double)> defined in <functional>. For reasons of efficiency, similar to the vector type, it is better to use a reference to the function rather than copying it, so we should use a reference, std::function<double(double)>& f when f is the input to a function. Unfortunately this suggests to C++ that we might be considering changing the function, which of course we are not allowed to do! So we have to declare the function argument constant, with const std::function<double(double)>& f to "promise" that we won't go changing f. We saw similar details when we looked at functions and call-by-reference earlier.

Of course, when we are using namespace std; we can eliminate std:: throughout.

2. Example: differentiation

// approxdiff.cpp (click here to download)

// Illustrates functional package and passing functions to as arguments
// This program contains a function that estimates derivatives
#include <iostream>
#include <functional>

using namespace std;

double derivative(const function<double(double)>& f, double x0, double h) {
  double h2 = h / 2;
  double lo = x0 - h2;
  double hi = x0 + h2;
  return (f(hi) - f(lo)) / h;
}

double g(double x) {
  return x * x / 2.0;
}

int main() {
  double x = 1;
  cout << "d/dx(x^2/2) at x=" << x << " is " << derivative(g, x, 1e-5) << endl;
  return 0;
}

Remark.

In fact the other variables can be declared const too: double derivative(const function<double(double)>& f, const double x0, const double h). Declaring all variables that are not varied as const is A Good Thing as it allows the compiler to spot programmer errors and to optimise code.

Remark.

The code above attempts to differentiate. How well it works depends on the small value of value h supplied by the user. Clearly h should be small to get reasonable accuracy. But there are other problems with small h as you will see. The actual code here makes sure that the derivative is not "biased" to one side or the other, and this is why $(f (x + h / 2) - f (x - h / 2)) / h$ is a better formula than the "obvious" $(f (x + h) - f (x)) / h$ . However, neither of these formulas are particularly good in many cases.

The value $h$ must be small, but when $h$ is small both formulas suffer significantly from catestrophic cancellation. Although differentiation is comparitively straightforward analytically (i.e. compared to integration), differentiation is particularly difficult to achieve numerically. On the other hand, it turns out that there are a great number of good techniques for numerical integration.

Of course, a small step size $h$ is necessary to give a good approximation to the derivative. Here is an example that shows how, contrary to popular belief, the value for $h$ can be so small that the approximation to the derivative is actually less accurate.

// diffsin.cpp (click here to download)

// Illustrates the numerical calculation of d/dx(sin x) at x = 1
// using different values of step size h
// uses <cmath> for the sin and cos functions

#include <iostream>
#include <functional>
#include <cmath>

using namespace std;

double derivative(const function<double(double)>& f, double x0, double h) {
  double h2 = h / 2;
  double lo = x0 - h2;
  double hi = x0 + h2;
  return (f(hi) - f(lo)) / h;
}

double g(double x) {
  return sin(x);
}

int main() {
  double x = 1;
  double h = 0.1;
  double ans = cos(x); // correct analytic solution for comparison
  double err = 0;
  cout.precision(10);

  cout << "exact value for d/dx(sin x) at x==1 is " << ans << endl;
  int n = 0;
  while (n < 30) {
    double app = derivative(g, x, h);
    double err = ans - app;
    if (err<0) { err = -err; }
    cout << "with h==" << h << " estimate is " << app << " error==" << err << endl;
    h = h/4.0; // make h go down by factor of 4
    n = n+1;
  }
  
  return 0;
}

Output of diffsin.cpp

Here is the program diffsin.cpp when run.

exact value for d/dx(sin x) at x==1 is 0.5403023059
with h==0.1 estimate is 0.540077208 error==0.0002250978217
with h==0.025 estimate is 0.5402882356 error==1.407026262e-05
with h==0.00625 estimate is 0.5403014265 error==8.793978465e-07
with h==0.0015625 estimate is 0.5403022509 error==5.496242883e-08
with h==0.000390625 estimate is 0.5403023024 error==3.435006501e-09
with h==9.765625e-05 estimate is 0.5403023057 error==2.14826823e-10
with h==2.44140625e-05 estimate is 0.5403023059 error==1.01905151e-11
with h==6.103515625e-06 estimate is 0.5403023059 error==1.473798861e-11
with h==1.525878906e-06 estimate is 0.5403023059 error==2.164179946e-11
with h==3.814697266e-07 estimate is 0.5403023056 error==2.693965051e-10
with h==9.536743164e-08 estimate is 0.5403023062 error==3.12680104e-10
with h==2.384185791e-08 estimate is 0.5403023073 error==1.476833322e-09
with h==5.960464478e-09 estimate is 0.5403023213 error==1.544667194e-08
with h==1.490116119e-09 estimate is 0.5403022468 error==5.905913403e-08
with h==3.725290298e-10 estimate is 0.5403023958 error==8.995247791e-08
with h==9.313225746e-11 estimate is 0.5403017998 error==5.060939698e-07
with h==2.328306437e-11 estimate is 0.540304184 error==1.878091821e-06
with h==5.820766091e-12 estimate is 0.5402946472 error==7.658651343e-06
with h==1.455191523e-12 estimate is 0.5403137207 error==1.141483499e-05
with h==3.637978807e-13 estimate is 0.5401611328 error==0.0001411730556
with h==9.094947018e-14 estimate is 0.5419921875 error==0.001689881632
with h==2.273736754e-14 estimate is 0.537109375 error==0.003192930868
with h==5.684341886e-15 estimate is 0.546875 error==0.006572694132
with h==1.421085472e-15 estimate is 0.46875 error==0.07155230587
with h==3.552713679e-16 estimate is 0.625 error==0.08469769413
with h==8.881784197e-17 estimate is 0 error==0.5403023059
with h==2.220446049e-17 estimate is 0 error==0.5403023059
with h==5.551115123e-18 estimate is 0 error==0.5403023059
with h==1.387778781e-18 estimate is 0 error==0.5403023059
with h==3.469446952e-19 estimate is 0 error==0.5403023059

Another reason why the formula $(f (x + h / 2) - f (x - h / 2)) / h$ is better than $(f (x + h) - f (x)) / h$ is seen using polynomial interpolation. Given the three data points $(x - h / 2, f (x - h / 2))$ , $(x, f (x))$ , $(x + h / 2, f (x + h / 2))$ we can interpolate a polynomial through these points using Lagrange's formula. Since there are three data points the polynomial that we fit will have degree at most 2. Its formula can be obtained and differentiated analytically, and the derivative evaluated at $x$ . The value one gets turns out to be exactly $(f (x + h / 2) - f (x - h / 2)) / h$ . (Try it, it is a page of messy but easy algebra!) So the "better" formula works exactly for all quadratics as well as all straight lines, whereas the other formula works exactly only for straight lines.

The previous remark shows why our method for differentiation was so good for $x^{2} / 2$ but poor for $\sin x$ . The reason is that, by Lagrange interpolation, the method is exact for quadratics but not exact for a function (such as $\sin x$ ) that is not quadratic. Of course $x^{2} / 2$ is quadratic, so the formula will be "perfect" for this one!

Remark.

Given all this, it seems to be a good idea to approximate the function $f$ to be differentiated at $x_{0}$ by a higher degree polynomial. The idea would be to take more data points, and use this to avoid having to go to very very small step sizes. Indeed this is often a very good idea, but doesn't avoid all the possible problems. So for example taking five equally spaced data points about $x_{0}$ ,

(x_{0} - 2 h, f (x_{0} - 2 h)), (x_{0} - h, f (x_{0} - h)), (x_{0}, f (x_{0})), (x_{0} + h, f (x_{0} + h)), (x_{0} + 2 h, f (x_{0} + 2 h))

we can fit these points to a quartic using Lagrange's formula and calculate the derivative of this quartic at $x_{0}$ . When we do so we get the following formula.

f' (x_{0}) \approx \frac{1}{h} (\frac{f (x_{0} - 2 h)}{12} - \frac{2 f (x_{0} - h)}{3} + \frac{2 f (x_{0} + h)}{3} - \frac{f (x_{0} + 2 h)}{12})

This formula gives, in many cases, a more accurate estimate of the derivative. One advantage is that, assuming that the function can be accurately approximated by a quartic, we can use a larger step size and hence reduce the effects of catastrophic cancellation. However, formulas like this do not work for all such functions $f$ .

// quarticdiff.cpp (click here to download)

// Illustrates functional package and passing functions to as arguments
// This program contains a function that estimates derivatives
#include <iostream>
#include <functional>
#include <cmath>

using namespace std;

/* estimate for derivative using quadratic approximation to f */
double derivative2(const function<double(double)>& f, double x0, double h) {
  double lo = x0 - h;
  double hi = x0 + h;
  return 0.5 * (f(hi) - f(lo)) / h;
}

/* estimate for derivative using quartic approximation to f */
double derivative4(const function<double(double)>& f, double x0, double h) {
  double xm2 = x0 - 2*h;
  double xm1 = x0 - h;
  double xp1 = x0 + h;
  double xp2 = x0 + 2*h;
  double h12 = 12.0*h;    // 12h/1
  double h23 = 3.0*h*0.5; // 3h/2
  return (f(xm2) - f(xp2)) / h12 + (f(xp1)-f(xm1)) / h23;
}

double g(double x) {
  return x * exp(x);
}

double gdiff(double x) {
  return (1.0 + x)*exp(x);
}

int main() {
  double x = 1; // the place to test
  double h = 0.5; // starting value of h;
  for (int n = 0; n<30; n++) {
    double d2 = derivative2(g,x,h);
    double d4 = derivative4(g,x,h);
    double d = gdiff(x);
    cout << "h==" << h << " \texact g'=" << d << " \terr_d2==" << abs(d2-d) << " \terr_d4==" << abs(d4-d) << endl;
    h = h/2.0;
  }
  return 0;
}

3. Summary

You have seen how in C++ 2011 it is possible to have functions as inputs to other functions. This is a useful technique that reduces repetition in coding.
You have seen how reducing step size in certain problems does - up to a point - increase accuracy. However reducing step size beyond a certain natural limit (determined by the problem in hand and the maximum accuracy available for double variables) actually decreases accuracy, and this effect is typical and pervasive.
You have seen how problems of catestrophic cancellation mean that numerical differentiation is particularly difficult to do numerically.

(It turns out that, in contrast to differentiation, numerical integration is much more successful, and there are a bewildering number of excellent methods for integration on a computer!)