Differentiation

Often one has to compute derivative of a function while doing numerical computation. If the form of the function is simple enough ( say sine function, exponential etc ) one knows the derivative analytically. Then the computation of the derivative is evaluation of the derivative function. That is straight forward. Many times the function form is quite complicated. In that case, the computation of analytical derivative could be involved. The expression of the derivative may be long and it may not be convinient to code the expression of the derivative. In that case one resorts to numerical computation of the derivative. Some cases, the function may not be given in a closed form and is evaluated numerically. In that case computation of analytical derivative is not possible. Finally, you may have a function given as a table of independent variable and the value of the function. Here also one has to resort to numerical differentiation.

The technique of numerical differentiation is based on the Taylor series expansion of a function:

$$f(x) = f(x_0) + ( x - x_0) \left . \frac{d f(x)}{dx} \right \vert... ...{(x - x_0)^2}{ 2 !} \left . \frac{d^2f(x)}{dx^2} \right \vert _{x = x_0} \cdots$$ (1)

If $x-x_0$ is small enough, one can neglect the higher order terms and define the approximate value of the derivative to be

$$\left . \frac{d f(x)}{dx}\right \vert _{x = x_0} = \frac{ f(x) - f(x_0) } { x - x_0 }$$ (2)

In fact, the derivative is defined as a limit of the expression above in the limit of $x$ approaching $x_0$. So, while computing the derivative on the computer, one has to evaluate the function f at two values of the argument close enough. But $x-x_0$ cannot be made arbitrarily small because of the round off errors.

Note that the error in the value of the derivative depends on $(x-x_0)^2$ and the second derivative of the function at $x_0$. If one can eliminate the second derivative somehow, one will be able to compute the derivative more accurately even for relatively large value of $x-x_0$. This can be done if we use the Taylor series for two values of x; $x_1 = x_0 + \Delta x$ and $x_2 = x_0 - \Delta x$. I leave it to you to show that the derivative then can be defined as

$$\left . \frac{d f(x)}{dx}\right \vert _{x = x_0} = \frac{ f(x_1) - f(x_2) } { x_1 - x_2 }$$ (3)

and then the error in the derivative is proportional to $(x_1 - x_2)^3 \left . \frac{ d^3 f(x)}{d x^3} \right \vert$.

Problem 1: Using these two definitions of the derivative, write a program to compute the derivative of known functions ( sine, cosine, exponential, logarithm, Bessel functions, power series etc ). Show that the second definition is superior.

One can extend the method of derivative calculation to high order derivatives. But it must be noted that the accuracy of the numerical values of the higher derivatives degrade very fast. This is because every time one computes a derivative there is loss of accuracy due to truncation. For higher order derivatives the error accumulates fast and soon one starts getting into noise.

Problem 2: Calculate first, second and third derivative of exponential numerically using one of the definitions given above. Show that the accuracy is least for third order derivative.

Note: While computing third order derivative, you need second order derivative at two values of the argument, this in turn means first order derivative is needed at four ( three if you are careful ) values of the argument and the function at eight ( or five ) values of the argument.

Many times one has data given as a function of an argument at fixed value. One then has to calculate the derivative. There are three different methods of computing the derivatives: methods using forward and backward differences and one using medians. The median method is the best since the values of the function used are close to the argument at which one is computing the derivatives. Further, if the values are given at equally spaced points, one can use the data such that even or odd derivatives are eliminated.