how to evaluate derivative of function in matlab?

Question

This should be very simple. I have a function f(x), and I want to evaluate f'(x) for a given x in MATLAB.

All my searches have come up with symbolic math, which is not what I need, I need numerical differentiation.

E.g. if I define: fx = inline('x.^2')

I want to find say f'(3), which would be 6, I don't want to find 2x

Answer 1

If your function is known to be twice differentiable, use

f'(x) = (f(x + h) - f(x - h)) / 2h

which is second order accurate in h. If it is only once differentiable, use

f'(x) = (f(x + h) - f(x)) / h     (*)

which is first order in h.

This is theory. In practice, things are quite tricky. I'll take the second formula (first order) as the analysis is simpler. Do the second order one as an exercise.

The very first observation is that you must make sure that (x + h) - x = h, otherwise you get huge errors. Indeed, f(x + h) and f(x) are close to each other (say 2.0456 and 2.0467), and when you substract them, you lose a lot of significant figures (here it is 0.0011, which has 3 significant figures less than x). So any error on h is likely to have a huge impact on the result.

So, first step, fix a candidate h (I'll show you in a minute how to chose it), and take as h for your computation the quantity h' = (x + h) - x. If you are using a language like C, you must take care to define h or x as volatile for that computation not to be optimized away.

Next, the choice of h. The error in (*) has two parts: the truncation error and the roundoff error. The truncation error is because the formula is not exact:

(f(x + h) - f(x)) / h = f'(x) + e1(h)

where e1(h) = h / 2 * sup_{x in [0,h]} |f''(x)|.

The roundoff error comes from the fact that f(x + h) and f(x) are close to each other. It can be estimated roughly as

e2(h) ~ epsilon_f |f(x) / h|

where epsilon_f is the relative precision in the computation of f(x) (or f(x + h), which is close). This has to be assessed from your problem. For simple functions, epsilon_f can be taken as the machine epsilon. For more complicated ones, it can be worse than that by orders of magnitude.

So you want h which minimizes e1(h) + e2(h). Plugging everything together and optimizing in h yields

h ~ sqrt(2 * epsilon_f * f / f'')

which has to be estimated from your function. You can take rough estimates. When in doubt, take h ~ sqrt(epsilon) where epsilon = machine accuracy. For the optimal choice of h, the relative accuracy to which the derivative is known is sqrt(epsilon_f), ie. half the significant figures are correct.

In short: too small a h => roundoff error, too large a h => truncation error.

For the second order formula, same computation yields

h ~ (6 * epsilon_f / f''')^(1/3)

and a fractional accuracy of (epsilon_f)^(2/3) for the derivative (which is typically one or two significant figures better than the first order formula, assuming double precision).

If this is too imprecise, feel free to ask for more methods, there are a lot of tricks to get better accuracy. Richardson extrapolation is a good start for smooth functions. But those methods typically compute f quite a few times, this may or not be what you want if your function is complex.

If you are going to use numerical derivatives a lot of times at different points, it becomes interesting to construct a Chebyshev approximation.