Compute a derivative using discrete methods

Question

I am looking for a method to compute a derivative using a discrete and fast method. Since now I do not know the type of equation I have, I am looking for discrete methods analog to the ones that we can find for the integral such as, the Euler method.

Answer 1

Automatic differentiation is the most accurate and conceptually awesome way of doing this kind of thing. Just a bit more complicated.

Answer 2

To do numerical differentiation, which is always an approximation, there are two common scenarios:

1. You have a way (algorithm, equation) to calculate the value of f(x) at any given x, or
2. You have the value of f(x) at a set of equally spaced values of x (f(1), f(1.5), f(2) etc):

1. If you have the algorithm, then you are best looking at Andrea Ambu's answer*:
   1. Start with the values of f(a+h) and f(a-h) and do
      f'(a) = { f(a+h)-f(a-h) } / 2h
      This is known as a central difference.
   2. Reduce the size of the offset h until f'(a) stops changing, as per Andrea's answer.
   • Be careful that you don't divide by 2h without checking it is larger than 0 or you'll get divide by zero errors.
2. If you have the value of the function at a set of values of f(x) = f(x_n) = f_n, then we can do something similar to the above method, but our accuracy will be limited by the values of x_n and the gaps between them.
3. The first approximation is just to use the same central difference operator as above;
   f'_n = (f_n+1 - f_n-1)/2h,
   where h is the equal spacing between the values of x_n.
4. The next is to use a five-point stencil, although it has only 4 coefficients as detailed on that wikipedia page. This uses the values from f_n-2 to f_n+2: f'_n = (-f_n+2 + 8f_n+1 - 8f_n-1 +f_n-2)/12h.
5. there are higher order approximations using more points, but they get increasingly difficult to calculate for diminishing returns, and become more unstable numerically.
6. **Note**: These finite-difference formulae rely on f being approximately the same shape as a polynomial in the range x_n-1≤ x ≤_n+1. For functions like sine waves, they can be quite bad at calculating the derivative to a good degree of accuracy.

* Andrea's answer uses the forward difference operator {f(a+h) - f(a)}/h instead of the central difference operator {f(a+h) - f(a-h)}/2h, but the forward difference operator is less accurate in numerical solutions.

Answer 3

I think you are looking for the derivative calculated in a point. If this is the case, here there is a simple way to do that. You need to know the derivative in a point, say a. It is given by the limit of the difference quotient for h->0:

You actually need to implement the limit function. So you:

• Define an epsilon, set it more small to be more precise, bigger to be faster
• calculate the difference quotient in a starting h, suppose h=0.01, store it in f1

• Now in a DO-WHILE loop:

  1- divide h by 2 (or by 10, the important thing is to make it smaller)
  2- calculate again the difference quotient with the new value of h, store this in f2
  3- set diff = abs(f2-f1)
  4- assign f1 = f2
  5- repeat from point 1 while (diff>epsilon)

• You can finally return f1 (or f2) as the value of your f'(a)

Remember: You are assuming the function is differentiable in a. Every result you will get will be wrong due of errors of the finite decimal digit your computer can handle, there is no escape from this.

Example in python:

def derive(f, a, h=0.01, epsilon = 1e-7):
    f1 = (f(a+h)-f(a))/h
    while True: # DO-WHILE
        h /= 2.
        f2 = (f(a+h)-f(a))/h
        diff = abs(f2-f1)
        f1 = f2
        if diff<epsilon: break
    return f2

print "derivatives in x=0"
print "x^2: \t\t %.6f" % derive(lambda x: x**2,0)
print "x:\t\t %.6f" % derive(lambda x: x,0)
print "(x-1)^2:\t %.6f" % derive(lambda x: (x-1)**2,0)

print "\n\nReal values:"
print derive(lambda x: x**2,0)
print derive(lambda x: x,0)
print derive(lambda x: (x-1)**2,0)

Output:

derivatives in x=0
x^2:         0.000000
x:       1.000000
(x-1)^2:     -2.000000


Real values:
7.62939453125e-08
1.0
-1.99999992328

The first time I've got "precise" value" because of using only the first 6 digits of the result, note I used 1e-7 as epsilon. The REAL calculated values are printed after that, and they are obviously mathematically wrong. The choose of how small epsilon is depends on how precise you want your results to be.

Answer 4

I hope this may be helpful NUMERICAL DIFFERENTIATION.

Answer 5

Formally, no. Either you are describing the (partial) derivitives of discrete functions or you are asking for a numerical method to approximate the (partial) derivatives of continuous functions.

Discrete functions don't have derivatives. If you review the epsilon-delta definition of a derivative, you will see that you would need to be able to evaluate the function close to the point you want the derivative at. That doesn't make sense if the function only has values at integer values of x, y and z. So there is no way to find the derivative of a discrete function for any value of fast.

If you want a numerical method exactly calculate the derivatives of a continuous function, you're out of luck as well. Numerical methods for derivatives are heuristic, not algorithmic. There is no numerical method which guarantees an exact solution. Fortunately, there exist many good heuristics. Mathematica uses a specialized version of Brent's principle axis method by default. I would recommend you use the GNU Scientific Library, which has a very good implementation of Brent's method. I owe my entire grade in one of my math courses to the GSL. The ruby bindings are pretty good if that's your thing. If necessary, most numerical differentiation libraries have a handful of different methods available.

If you really want, I can whip out some sample code. Let me know.

Answer 6

If the function is linear as you have indicated then the derivatives are trivial. The derivative with respect to 'x' is 'a'; the derivative with respect to 'y' is 'b' and the derivative with respect to 'z' is 'c'. If the equation is of a more complex form and you need a formula representing the solution rather than an empirical solution, then please submit the more complex form of the equation.

Regards

Answer 7

There is quite a bit of theory (and established practice) in calculating numerical ("finite") derivatives. Getting all the details correct, such that you believe the result, is not trivial. If there's any way you can get the analytical derivative of the function (using pen and paper, or a computer algebra system such as Maple, Mathematica, Sage, or SymPy), this is by far the best option.

If you can't get the analytical form, or you don't know the function (just it's output), then numerical estimation are your only option. This chapter in Numerical Recipies in C is a good start.

Answer 8

Short of using a symbolic math language, like Maple, the best you can do is to approximate the derivative at various points. (And then interpolate, if you need a function.)

If you've already got the function you want to use, then you should use the backward divided-difference forumla, and Richardson extrapolation to improve your error.

Also keep in mind that these methods work on functions of one variable. However, partial derivatives of each variable treat the other variables as constants.

Answer 9

A simple method is to compute the change in f over a small value for each point of the derivative you're interested in. For example, to compute ∂f/∂x, you could use this:

epsilon = 1e-8
∂f/∂x(x, y, z) = (f(x+epsilon,y,z) - f(x-epsilon, y, z))/(epsilon * 2);

The other partials would be similar in y and z.

The value chosen for epsilon depends on the contents of f, the precision required, the floating point type used, and probably other things. I suggest you experiment with values for it with functions you're interested in.

Answer 10

I'll assume that your function is more complex than the simple one you posted, because the closed-form solution is far too simple.

When you use the word "discrete" it makes me think you need "finite differences". You'll need some discretization to calculate the approximation.

Df/Dx ~ (f2-f1)/(x2-x1), etc.