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.

  • 1
    $\begingroup$ The wikipedia entry Numerical differentiation contains a lot of important information and several methods that are usually implemented. $\endgroup$
    – AccidentalFourierTransform
    Dec 14 '16 at 17:56
  • $\begingroup$ Hi. The Wikipedia entry talks mostly about the finite difference method and the related problem with it regarding the rounding off error of the machine. That is the kind of problem you run into if your step size is too small. I also ran into the same problem in my original calculation (actually another version of it). There is an "ideal" step size for which the error is minimum and going either way increases the error. There some text of quadrature method so may be someone can point into that direction. $\endgroup$
    – Shaz
    Dec 14 '16 at 18:02
  • $\begingroup$ Doesn't this question belong on a different site, like Mathematics? $\endgroup$
    – march
    Dec 14 '16 at 18:39
  • $\begingroup$ There is something wrong with your implementation. These errors are far too large to be explained by a problem with the method. $\endgroup$ Dec 14 '16 at 18:39
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Dec 16 '16 at 9:26

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.


  • 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.


  • 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.


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


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


  • 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).
  • $\begingroup$ Hi, thanks for the comment. The problem which appears that there is no universal $h$ which is adequate for all points on the grid. So, to say, if I find a certain $h$ which gives a good enough answer for a small domain of $f(s)$ will not work for another domain, one which is much closer to $x=0$. $\endgroup$
    – Shaz
    Dec 15 '16 at 5:33
  • $\begingroup$ That's why I wrote "adaptively", meaning you chose (automatically) $h$ differently in different parts of the grid. Eventually you will run into other problems (rounding, precision) but this can't be helped if you are trying to do numerics with ugly functions. $\endgroup$
    – user1583209
    Dec 15 '16 at 7:07
  • $\begingroup$ I am beginning to realize this as I delve in the problem further with other functions that I encounter in my original calculation. I am trying to write a code to calculate the Ricci tensor and Ricci scalar for an arbitrary given metric tensor. I cant imagine that I am the only person that has encountered this problem. The situation is that apart from Minkowski metric in cartesian coordinate, all other matrix functions are non-linear (involving inverse of sine and radius). My code gives me the correct answer for cartesian coordinate but that a very trivial answer. $\endgroup$
    – Shaz
    Dec 15 '16 at 7:33
  • 1
    $\begingroup$ Your problem is not non-linearity, but the divergence. $\endgroup$
    – user1583209
    Dec 15 '16 at 9:21
  • 2
    $\begingroup$ @Shaz: In my understanding, "non-linear" means that it is not a linear function (i.e. not a straight line if you plot it). What you have here is a singularity, an infinite discontinuity, which is something different. $\endgroup$
    – user1583209
    Dec 16 '16 at 9:33

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy