I'm trying to implement a finite difference derivative algorithm, as I'm learning how to perform numerical analysis on computational physics. Now I'm following the book found here. I'm struggling to implement the algorithm used. Here's a bit of background:

As you may know, the function is divided into finite segments called $h$. And each discrete term of the function is labelled as $x_i$. However, the book also defines an arbitrary value called $\zeta$ where $\zeta\in|x_i,x_{i+1})$. It also defines another variable $\varepsilon$ such that $\zeta=x_i+\varepsilon h$. Note that $\varepsilon\in[0,1)$. I know that they're in between the discrete values, but I can't tell what they are.

And to show some conventions used, $f^{(n)}(x_i)\equiv f^{(n)}_i$; and $f^{(n)}_{i+\varepsilon}=f^{(n)}(\zeta)$.

Now if you're familiar with finite difference differentiation, you may know that in order to perform it, you may expand the general term using a Taylor series to change the order of the truncation error. So what the book did was to expand it to a certain number of terms, and played around with the forward difference, backward difference, and central difference approximation terms to get the following equation (I omitted the details to keep the question short):

$f'_i=\frac{1}{6h}(f_{i-1}+8f_{i+\frac{1}{2}}-8f_{i-\frac{1}{2}}-f_{i+1}) + \frac{h⁴}{360}f^V_{i+\varepsilon}$

Now when I tried to punch it into my computer, the approximation I got was way off. I didn't enter second term that has the fourth power of $h$ as I can't tell what $\varepsilon$ or $\zeta$ is. Here is an example code of when I try to approximate $x ^ 2$ when $x=2$:

let derivative = 1.0/(6.0*h) * (1.0+8*(2.5*2.5)-8.0*(1.5*1.5)-9.0)
// The above line gives an approximation of 400, when it should be close to 4

What am I doing wrong here? If there's anything more you need to know, please say so. I'm sorry that the question is so long and wordy.


Can you be explicit about what you used for $h$? It looks like you used $h = 1$? That's likely going to give an inaccurate result. The quality of the approximation gets better as $h$ gets smaller. Try a sequence of decreasing values of $h$. If you take $h$ "too large" you still get a number out, but it's not a useful number.

EDIT After a couple of exchanges in the comments, it's clear that you're not using the same value of $h$ in the factor $1/(6h)$ in the front of the expression that you're using in the evaluation of the function. You claimed to use $h=0.01$. Looking at the values that you used $f_{i-1} = 1$ (as well as the other specific values that you showed for the other evaluations), however, shows that you were actually spacing the function evaluations by $h=1$. You need to use the same value of $h$ throughout.

Also, as I pointed out in the comments, the value of a good $h$ varies a lot from problem to problem. Knowing that some specific value works on a problem that you solved before does not necessarily mean that it will work on the next problem.

Edit 2 The variables $\zeta$ and $\epsilon$ are giving the error estimate based on the Lagrange form of the remainder for the Taylor series that you're using to construct the finite difference. In practice you will just drop that term when you do most calculations, but if you want to do an analysis of numerical errors introduced you will want it. Even then you usually only use the factor in front that's proportional to $h$ since you usually want to know how the numerical errors scale with the discretization length. If you wanted a different type of bound, you might need the whole term. In general you don't know the exact values of $\zeta$ and $\epsilon$ and would trying to bound not calculate that term exactly.

See https://en.wikipedia.org/wiki/Finite_difference_method#Accuracy_and_order

  • $\begingroup$ It was 0.01. It was good enough for other approximations to a fine degree. It makes no sense for this to be equal to 400, when it should be 4 only. $\endgroup$ – user3894009 Sep 17 '15 at 19:19
  • $\begingroup$ To add to that, I'm still not 100% sure what $\zeta$ and $\varepsilon$ are used for. But as I said, when computing it using the normal definition of the derivative, but using $h=0.01$ instead of a limit, it worked very well. $\endgroup$ – user3894009 Sep 17 '15 at 19:23
  • 1
    $\begingroup$ @user3894009 You need to be consistent in your value of $h$ in both the factor $1/(6h)$ and the points where you evaluate $f$. In your example it looks like you put $f_{i-1}=1$. Since you are evaluating at $x=2$, that looks like you have $h=1$ in that part. Likewise it looks like you have $f_{i/2+1} = (2.5*2.5)$, which again seems to imply that you were using $h=1$. If you have $h=0.01$, then you should have, e.g., $f_{i+1}=(2.01)^2$. $\endgroup$ – Brick Sep 17 '15 at 19:27
  • $\begingroup$ @user3894009 In this case, I agree that $h=0.01$ should give a more accurate answer than you're getting (assuming that you do the calculation correctly), but also keep in mind that "good" values of $h$ are highly dependent on the problem. A value that worked on one problem maybe terrible on another. $\endgroup$ – Brick Sep 17 '15 at 19:30
  • $\begingroup$ Good catch. I was adding 1 instead of $h$. I didn't notice it. Thanks, it works now! $\endgroup$ – user3894009 Sep 17 '15 at 19:43

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.