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.
Answer
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
No comments:
Post a Comment