9.1 Numerical Differentiation
How can we find a good approximation to the derivative of a function?
The obvious approach is to pick a very small \(d\) and calculate \(\frac{f(x+d)-f(x)}{d}\), which looks like the definition of the derivative. Actually, this is not a great idea.
Why?
The problem is that your means of calculating are not infinitely accurate, particularly if \(f(x)\) or \(f(x+d)\) are irrational numbers. This implies that there will sometimes be small errors in your evaluations. \(f(x+d)\) and \(f(x)\) will differ from each other by something like \(f'(x)d\) when \(d\) is small, but your calculation errors in them will be roughly independent of \(d\). As a result, as you make \(d\) smaller, the ratio of your error to \(f'(x)d\) grows. Dividing the result by a very small \(d\) is the same as multiplying it by a very large \(\frac{1}{d}\), and that amplifies the error. When \(d\) becomes smaller than the size of your calculational errors, the estimate you get for the derivative will be mostly calculational error and will tell you very little about \(f'(x)\). The kind of error that arises in this way is usually called round-off error.
The upshot of this is that you really want to compute the derivative using only relatively large values of \(d\).
Is this possible?
The answer is yes! And doing it is lots of fun.
How?
Here is the basic idea: Suppose your function \(f\) is not only differentiable, but its derivative is also differentiable at argument \(x\), and so is its derivative, etc.
If so, the value of \(f(x+d)\) can be described by a power series,
\[f(x+d) = f(x) + f'(x)d + f''(x)\frac{d^2}{2} + f^{(3)}(x)\frac{d^3}{3!} + ...\]
(here \(f^{(j)}(x_0)\) means the \(j^{th}\) derivative of \(f(x)\) evaluated at \(x = x_0\).) To prove this, calculate the derivatives of both sides with respect to \(d\) at \(d = 0\), and the second derivatives etc.)
We want \(f'(x)\), so we want to get rid of the \(d^2\) and further terms on the right.
If we form \(\frac{f(x+d)-f(x)}{d}\) we will get \(f'(x) + f''(x)\frac{d}{2} + f^{(3)}(x)\frac{d^2}{3!} + ...\) and there is an error term of \(f''(x)\frac{d}{2}\) that is proportional to \(d\).
On the other hand, if we instead form \(\frac{f(x+d)-f(x-d)}{2d}\) all the terms in the series above that have even powers of \(d\) disappear, and we get \(f'(x) + f^{(3)}(x)\frac{d^2}{3!} + f^{(5)}(x)\frac{d^4}{5!} + ...\) The error term in this expression is proportional to \(d^2\).
This is already a big improvement over the obvious estimate. Here the error decreases as \(d^2\) rather than as \(d\) as you reduce \(d\). The wonderful thing is, we can do even better, by eliminating the \(f^{(3)}\) term, and then the \(f^{(5)}\) term, and so on, as far as we want to go.
How can we do that?
Well, you can combine your evaluation of \(\frac{f(x+d)-f(x-d)}{2d}\) with one of \(\frac{f(x+2d)-f(x-2d)}{4d}\). The second of these will have the same \(f'(x)\) term in the power series, but \(4\) times more of the \(d^2\) term. Thus if we form \(4\) times the first of these and subtract the second, we will end up with three times \(f'(x)\), no \(d^2\) term at all, and only correction terms of order \(d^4\) and higher.
Thus if we define \(C(d)\) to be \(\frac{f(x+d)-f(x-d)}{2d}\), the combination \(\frac{4C(d)-C(2d)}{3}\) will yield \(f'(x)\) plus an error coming from the fifth derivative term in the original expansion rather than the third, and that error term will be proportional to d4 in our computation (plus terms proportional to \(d^6\), \(d^8\), etc.).
Call this combination \(D(d)\); then similarly, \(\frac{16D(d)-D(2d)}{15}\) (call this \(E(d)\)) will yield \(f'(x)\) plus an error that comes from the seventh derivative, and is proportional to \(d^6\). And you could go on to form \(256\) times this \(E(d)\) minus its value for twice \(d\) all divided by \(255\) to get an expression whose error will be proportional to \(d^8\).
What this means is that dividing \(d\) by \(2\) will reduce the error in this last estimate \(E(d)\) by a factor of \(2^8\) which is \(256\).
This looks like a mess.
But it isn't. It is very easy to do all this on a spreadsheet, and you can see what happens to each of the estimates above as you reduce \(d\) by a factor of \(2\) successively, for any function you can write down, and any argument \(x\).
Not only that, you can change the argument by changing only one entry, and change the function by changing only one entry and doing some copying.
OK, how?
We will set this up using the function \(\sin(x)\) to be specific
Preliminaries:
1. Put Calculating \(f'(x)\) in A1
2. Put the name of your function in A2
3. Put dstart in A3 and your starting value for \(d\) in B3 (I put \(1\) in B3)
4. Put the letter \(x\) in A4 and your value of \(x\) in B4 (I put \(1\) there also)
5. Label columns in row 5 as follows: in A5 put \(d\), in B5 \(x\), in C5 \(x+d\), in D5 \(x-d\), in E5 \(f(x)\), F5 \(f(x+d)\), G5 \(f(x-d)\), H5 \(C(d)\), in I5 \(D(d)\), in J5 \(E(d)\)
Setup
Now in A6 enter =B3, in B6 enter =B$4, in C5 enter =A6+B6, in D6 =B6-A6, in E6 enter =f(B6). For example =tan(B6)
Then copy E6 to F6 and G6. In H6 enter =(A6-G6)/2/A6
Copy B6 through H6 down their columns to row 50
Now enter in A7 =A6/2 and copy A7 down to row 50. In I7 enter =(4*H7-H6)/3 and copy I7 down to row 50
Finally in J8 enter =(16*I8-I7)/15 and copy J8 down to row 50.
What is all this?
Column A will contain the difference d used in your computation. It will divide by \(2\) in going from one line to the next. Column B contains the value of \(x\) at which you are finding the derivative, C and D contain \(x+d\) and \(x-d\). E contains \(f(x)\) and F and G contain \(f(x+d)\) and \(f(x-d)\). H contains the estimate \(\frac{f(x+d)-f(x-d)}{2d}\), column I contains the improvement obtainable by taking four times the H estimate, subtracting the similar estimate for \(d\) replaced by \(2d\), and dividing this difference by \(3\). J contains the improvement obtained by similarly subtracting the \(2d\) estimate in I from \(16\) times the \(d\) estimate there, and dividing by \(15\).
Here is the result for the function \(\tan(x)\) at \(x = 1\).
Notice that column E is accurate to \(9\) places when \(d\) is around \(\frac{1}{100}\).
To change \(x\) you need only put the value you want in B4. To change functions you enter the new function in E6 with variable B6, copy to F6 and G6, and copy columns E, F, and G down to the \(50^{th}\) row. The other columns need not be changed at all.
Exercises.
9.1 Set this up yourself. What is the value of d when the difference between E and F reaches the ten digit accuracy?
9.2 Try finding \(\tan(1)\) using \(\frac{f(x+d)- f(x)}{d}\) instead of E as above. For what \(d\) value do you reach the correct answer to ten digit accuracy?
9.3 Find a function and value for which H above does not get the answer to ten digit accuracy.