Numerical differentiation of non-linear functions

Question

In my efforts to write a code for a calculation I have encountered a problem of numerically differentiating a non-linear function at different points on a grid. I used the simple forward finite difference method to differentiate the function but the errors were huge in case when I know the analytical solution of the calculation.

So to better understand the problem I numerically differentiated $f(x)=1/x$ since the analytical solution is known to be $f'(x)/=-1/x^2$. So I tried to calculate the error between numerical and analytical derivative of $f(x)$ and plot the error against the step size for various values of $x$. The result is attached as a picture(Both axes are logarithmic!).

It appears from the graph the error is significant event for the smallest value of stepsize when $x$ is small. That makes sense to me now because a straight line approximation would fail as the function becomes more and more non-linear with decreasing $x$.

So my question is, is there a better way of numerically differentiating a non-linear function on the domain of the function other than the finite difference method? The function(s) that I need to numerically differentiate are multivariable functions but the method can be easily extended from univariable functions to multivariable function I suppose.

Answer 1

You are missing a minus sign in the derivative, $f'(x)=-\frac{1}{x^2}$.

The error to first order in $h$ is

$$\frac{f(x+h)-f(x)}{h} - f'(x)\approx \frac{h}{x^3}$$ which fits with your data as far as I can see. If $x$ is small, this can be a very big number for fixed $h$. So your problem is purely mathematical and not due to mistakes in your implementation or to numerical rounding problems.

The problem is of course that $f(x)$ is diverging as $x\to 0$.

You can deal with problems of this type by (suggestions are not limited to your specific case of $1/x$ or to numerical differentiation, but generally to dealing with "ugly" not well-behaved functions numerically):

• Transforming the function to a nicer looking function (without divergences, etc); In this specific case you could multiply with $x$ to get $g(x)=f(x)x$, differentiate $g(x)$ numerically and use that $f'(x)=(g'(x)-f(x))/x$. Another possible transform could be, that if you are dealing with a complicated function g(x) of which you know the asymptotic behavior at $x=0$ to be like $1/x$ (or any other simple function you can handle analytically), you could instead work with g(x) - 1/x which would not have the singularity.

or

• pick a sufficiently small $h$, perhaps adaptively (smaller as you get closer to the divergence); Eventually, if you are getting too close to 0, you will run into problems with precision or rounding, so this method alone will not get you all the way to the singularity, but if you are fine with values a bit away from the singularity it could suffice.

or

• If it is a real world problem (from physics,...), analyze it from that perspective. Usually things don't go to infinity in the real world, so perhaps the quantity you are dealing with is not the right quantity to analyze.

or

• If this is part of a function, perhaps another part cancels this singularity!?

or

• See if shifting the singularity with transforms $x\to x_0 - x$ or similar helps.

or

• See whether other numerical methods are more suited to your task. (E.g. in your case you could try the central difference instead of forward difference).

Answer 2

The answer to this question depends on what kind of functions you are actually interested in differentiating. That said, I will give a basic, general answer here.

As already pointed out, you are using a finite difference formula

$$f'(x) \approx \frac{f(x+h)-f(x)}{h}.$$

that is accurate to first order (i.e., the error is proportional to $h$). There are two simple ways to improve the accuracy. The first, as you know, is to decrease $h$. The only issue with this is that for very small $h$, roundoff error will become significant (you have not come close to such small values in your plot). The other general way to get more accurate approximations is to use a formula with a higher order of accuracy. For instance, the centered formula

$$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}$$

is accurate to second order (i.e., the error is proportional to $h^2$). For a better understanding of finite difference approximations of derivatives, and their accuracy, read any introductory numerical analysis book, or even just Wikipedia. I particularly like chapter one of LeVeque's text.