7.3. Numeric DifferentiationΒΆ

Computing a derivative numerically, is typically not difficult, but the accuracy of the calculation may be limited by the sampling rate of the data or by the limits of the computer to process very small numbers.

Calculus defines the derivative of y(x) with respect to x as

\frac{dy(x)}{dx} = \lim_{\Delta\,x \rightarrow 0}
\frac{\Delta\,y(x)}{\Delta\,x}.

For discretely sampled data, this becomes an approximation.

\frac{dy_i}{dx} \approx \frac{y_{i+1} - y_i}{h}

The variable h represents the distance with respect to x of the samples.

For sampled data values, the MATLAB function diff returns the difference between consecutive values, dy_dx = diff(y)/h;. Note the size of the array returned by diff is one element smaller than the array passed to it.

To take the derivative of a function, first sample the data at intervals of h; and then use diff as with pre-sampled data. The challenge is to determine an appropriate value for h. It should be small enough to give good granularity. The evaluation of small enough depends entirely on the data being evaluated. However, it should not be so small as to require excessive computation or to cause excessive rounding errors from the division by h. As rule of thumb, h should not be much smaller than 10^{-8}. Some testing of different values of h may be required for some applications.

Try this example code and see how changing the value of h, affects the accuracy of the numeric derivative.

% Derivative example, Compute deivative of cos(x)

% h = 0.1;            % sample step size, please experiment with h values
h = pi/40;
x = 0:h:2*pi;
y = cos(x);
dy_dx = diff(y)/h;
x1 = x(1:end-1);      % note one less to match dy_dx
actual_dy = -sin(x1);
error = abs(dy_dx - actual_dy);
me = max(error);
ex = x1(error == me);

p1 = plot(x1, actual_dy, 'k', x1, dy_dx, 'r');
title('Derivatives of cos(x)');
xticks(0:pi/2:2*pi);
xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});
line([ex ex], [-1 1]);
text(ex+0.1, -0.4, 'Max error location');
legend(p1,{'analytic derivative','numeric derivative', },...
    'Location', 'northwest');
disp(['Max error = ',num2str(me)]);