Numerical Derivatives

Forward Difference Derivative:

A simple approximation for this is to simply evaluate the above expression for a small, but finite, h.

 d         y(x+h) - y(x)
---- y =  ---------------
 dx              h

This is known as the forward difference derivative. Given n (x,y) points, we can then evaluate y', (or dy/dx), at n-1 points using the above formula.

The forward difference derivative can be turned into a backward difference derivative by using a negative value for h. Alternatively, many consider the two point formula as a method for computing not y'(x), but y'(x+h/2), however this is technically a three point derivative analysis.

Error Analysis:

Notice two things about this approach. First, we have approximated a limit by an evaluation which we hope is close to the limit. Second, we perform a difference in the numerator of the problem. This gives us two sources of error in the problem. Large values of h will lead to error due to our approximation, and small values of h will lead to round-off error in the calculation of the difference.

The error due to the finite difference approximation is proportional to the width h. We refer to the error as being of the order of h.

Three Point Formulas:

Three Point Methods allow us to use information from both before and after the point to be evaluated to better evaluate the derivative.

The most common three point method is an average of a forward and backward difference derivative. The approximation errors in the forward and backward difference schemes cancel, leaving approximation error of the order h², that is, the error is proportional to the grid width squared (remember, for h<1, h² is less than h).

               y(x+h) - y(x-h)
     y'(x) =  -----------------
                     2*h

Notice how this reduces to the interpretation of the two point formula more accurately representing the value of y'(x+h/2) if one substitutes h/2 for h in the above formula. Try this in the applet by switching back and forth from the two point to the three point first derivative calculation.

This method will allow you to solve for y' at n-2 points. The endpoints cannot use this formula, because we do not know y(x-h) for our first point, or y(x+h) for our last point.

However, a different three point method exists which allows us to evaluate the value of y' at the endpoints.

            -3*y(x) + 4*y(x+h) - y(x+2h)
     y'(x) = ------------------------------
                        2*h

Notice that this can be used for the final endpoint by replacing h with -h.

Error Analysis:

Both the standard and endpoint three point formulas have approximation error on the order of h².

Second Derivatives:

The second derivative of a function is simply the derivative of the derivative.

          d           
y''(x) = ---- y'(x) 
          dx         

or

          d2
y''(x) = ----    y(x)
          dx2

Forward Difference:

The simplest way to calculate this is to simply apply the forward difference formula at n points to get y'(x) at n-1 points, and then apply the same formula to y'(x) at n-1 points to get y''(x) at n-2 points.

Three Point Formula:

A three point formula can be constructed which uses the difference in results of the forward and backward two point difference schemes, and computes a three point derivative of that to get the second derivative.

          y(x-h) - 2*y(x) + y(x+h)
y''(x) = --------------------------
               h2

Endpoint Evaluation:

The above formula suffers from the same problem at endpoints as the three point formula for the first derivative. One approach would be to determine the first derivative using the endpoint methods above, and then to compute the second derivative from the first derivative in the same fashion.

Uneven Grid Spacing:

The above formulas all assume equal grid spacing. Equal grid spacing makes it easier to achieve higher degrees of precision in numerical derivative calculation, and should be used when possible.

The standard way of viewing this is by assuming a Taylor's series, which approximates y(x) in terms of the first, second, third, and subsequent derivatives. The idea is that in general, the rate of change of the function y(x) says more about what is happening at a specific instant than the rate of change of the rate of change, or the rate of change of the rate of change of the rate of change. As you move to higher order derivatives, the impact lessens.

This is written as:

y(x + h) = y(x) + h * y'(x) + 
           h2/2 y''(x) + 
           h3/6 y'''(x) + ... +
           hn/n! y(n)(x) 
           + ...

Two things should become apparent looking at this. First, as long as h<1, higher order terms will, in general, be small. Second, we can get y(x - h) by substituting -h for h in the expansion, yielding

y(x - h) = y(x) - h * y'(x) +
           h2/2 y''(x) - 
           h3/6 y'''(x) + ... +
           (-1)nhn/n! y(n)(x)
           + ...

Now, we mentioned above the results for even grid spacings in terms of "averages" of the forward and backward derivatives. Notice that by adding the above expressions for y(x + h) and y(x - h) that we get the exact same solution (plus error terms on the order of h⁴)!

However, if the grid spacing is not even, then we are no longer adding y(x + h) and y(x -h), but y(x + h) and y(x - g) where g is not equal to h. All of our errors which cancelled before no longer cancel out! Also, we need to know the first derivative in order to calculate the second derivative for unequal grid spacing. We get higher order error both in the fact that our order h³ error no longer cancels out (the y'''(x) term), but also in having to numerically calculate y'(x), which also no longer cancels out.

References:

Burden, R. L. and Faires, J. D., Numerical Analysis, Fifth Edition. pp. 156-167. PWS Publishing Co. Boston, MA, 1993.