Finite difference differentiation formula

Question

I'm trying to understand how the co-efficients of finite differences are calculated.

In particular I'm interested in the first derivative for a uniform grid of unit width.

I found this document which tabulates the co-efficients but seems a little sketchy on their derivation.

This paper mentions on page 24 appendix B, a formula for calculating the derivative, but I'm not sure about the derivation, or how to calculate C.

(I also have some doubts about the limits on the big sigma. Shouldn't they go from -(n-1)/2 to (n-1)/2 ?

I found a formula on the wikipedia page

$\delta^n_h[f](x)=\Sigma_{i=0}^{n} (-1)^iC_i^nf(x+(\frac{n}{2} -i)h)$

While this formula does give me the first derivative, it does not allow me to choose the number of points used in the calculation.

I'm specifically looking for the derivation of the formula for j=17 on page 937 of this document

Does anyone know what formula is used for this?

Does anyone know how this is derived?

Answer 1

The formulas have to be correct for the exponential function $f(x)=\exp(x)$. Conversely, by Laplace-transform or Taylor shift arguments, if the formula is true for the exponential, then it holds for all functions with derivatives of sufficiently high order. Thus $$ 1=\frac{a_1·(e^h-e^{-h})+...+a_k·(e^{kh}-e^{-kh})}{x}+O(h^p) $$ has to be true for all $h$. Dividing out the common factor $(e^h-e^{-h})=2·\sinh(h)=4·\sinh(h/2)·\cosh(x/2)$ results in $$ h=\sinh(h)·(b_0+b_1·(2\sinh(h/2))^2+b_2·(2\sinh(h/2))^4+…+b_{k-1}·(2\sinh(h/2))^{2(k-1)})+O(h^{p+1}) $$ and with $z=2\sinh(h/2)$, $h=2·\text{Arsinh}(z/2)$ $$ 2·\text{Arsinh}(z/2)=z·\cosh(\text{Arsinh}(z/2))·(b_0+b_1·z^2+b_2z^4+…+b_{k-1}·z^{2(k-1)})+O(z^{2k+1}) $$ one can read of the wanted coefficients by comparing truncated power series. $$ B_k(z)+O(z^{2k+1})=\frac{2·\text{Arsinh}(z/2)}{\cosh(\text{Arsinh}(z/2))}. $$ The relevant power series starts as

> as:=Arsinh(Z/2+O(Z^40)); bb:=2*as/Cosh(as); bb;
    Z - 1/6*Z^3 + 1/30*Z^5 - 1/140*Z^7 + 1/630*Z^9 - 1/2772*Z^11 
      + 1/12012*Z^13 - 1/51480*Z^15 + 1/218790*Z^17 - 1/923780*Z^19 
      + 1/3879876*Z^21 - 1/16224936*Z^23 + 1/67603900*Z^25 
      - 1/280816200*Z^27 + 1/1163381400*Z^29 - 1/4808643120*Z^31 
      + 1/19835652870*Z^33 - 1/81676217700*Z^35 + 1/335780006100*Z^37 
      - 1/1378465288200*Z^39 + O(Z^40)

There the coefficient of $Z^{2k+1}$ is also the coefficient of the error term before $h^{2k}·f^{(2k+1)}(x_v)$ for the formula using values $f_{-k},…,f_k$.

Undoing all the transformations and using $F=\exp(x/2)$, $\cosh(x/2)=(F+F^{-1})/2$, $\sinh(x/2)=(F-F^{-1})/2$ so that $F^{2j}$ stands for $f_j=f(x+j·h)$ one gets $$ a_1·(F^2-F^{-2})+…+a_k·(F^{2k}-F^{-2k})=(F+F^{-1})/2·B_k(F-F^{-1}) $$

> [ Evaluate(Truncate(bb+O(Z^(2*k))),(F-F^-1))*(F+F^(-1))/2 : k in [1..9] ];                    
[
    -1/2*F^-2 + 1/2*F^2,

    1/12*F^-4 - 2/3*F^-2 + 2/3*F^2 - 1/12*F^4,

    -1/60*F^-6 + 3/20*F^-4 - 3/4*F^-2 + 3/4*F^2 - 3/20*F^4 + 1/60*F^6,

    1/280*F^-8 - 4/105*F^-6 + 1/5*F^-4 - 4/5*F^-2 + 4/5*F^2 
        - 1/5*F^4 + 4/105*F^6 - 1/280*F^8,

    -1/1260*F^-10 + 5/504*F^-8 - 5/84*F^-6 + 5/21*F^-4 - 5/6*F^-2 
        + 5/6*F^2 - 5/21*F^4 + 5/84*F^6 - 5/504*F^8 + 1/1260*F^10,

    1/5544*F^-12 - 1/385*F^-10 + 1/56*F^-8 - 5/63*F^-6 + 15/56*F^-4 
        - 6/7*F^-2 + 6/7*F^2 - 15/56*F^4 + 5/63*F^6 - 1/56*F^8 
        + 1/385*F^10 - 1/5544*F^12,

    -1/24024*F^-14 + 7/10296*F^-12 - 7/1320*F^-10 + 7/264*F^-8 
        - 7/72*F^-6 + 7/24*F^-4 - 7/8*F^-2 + 7/8*F^2 - 7/24*F^4 
        + 7/72*F^6 - 7/264*F^8 + 7/1320*F^10 - 7/10296*F^12 
        + 1/24024*F^14,

    1/102960*F^-16 - 8/45045*F^-14 + 2/1287*F^-12 - 56/6435*F^-10 
        + 7/198*F^-8 - 56/495*F^-6 + 14/45*F^-4 - 8/9*F^-2 + 8/9*F^2 
        - 14/45*F^4 + 56/495*F^6 - 7/198*F^8 + 56/6435*F^10 
        - 2/1287*F^12 + 8/45045*F^14 - 1/102960*F^16,

    -1/437580*F^-18 + 9/194480*F^-16 - 9/20020*F^-14 + 2/715*F^-12 
        - 9/715*F^-10 + 63/1430*F^-8 - 7/55*F^-6 + 18/55*F^-4 
        - 9/10*F^-2 + 9/10*F^2 - 18/55*F^4 + 7/55*F^6 - 63/1430*F^8 
        + 9/715*F^10 - 2/715*F^12 + 9/20020*F^14 - 9/194480*F^16 
        + 1/437580*F^18
]

Complete Magma commands to get the output above, e.g. in the online calculator

P<Z>:=PowerSeriesRing(Rationals());  
L<F>:=LaurentSeriesRing(Rationals());
Arsinh := func<y|Log(y+Sqrt(1+y^2))>;  
as:=Arsinh(Z/2+O(Z^40)); bb:=2*as/Cosh(as); bb;    
apr := [ Evaluate(Truncate(bb+O(Z^(2*k))),F-F^-1)*(F+F^(-1))/2 : k in [1..19] ]; apr;

Answer 2

There is a Maple toolkit (FD) that computes the finite difference stencil expressions for given points and derivative upto arbitrary accuracy. Check out section 2.1 (computing the FDA expressions) for the details of how to compute the coefficients. The FD toolkit computes the coefficients for you.

There is a simple front-end function in FD:

Sten(diffexpr,[points])

which computes the stencil using the points, [points]. For example see the following demonstration of computing the forward and centered FDA for the first and second derivatives:

> Sten(diff(f(x),x),[0,1,2]);
                           -3 f(x) + 4 f(x + h) - f(x + 2 h)
                       1/2 ---------------------------------
                                           h
 > Sten(diff(f(x),x,x),[-1,0,1,2,3]);
             11 f(x - h) - 20 f(x) + 4 f(x + 2 h) + 6 f(x + h) - f(x + 3 h)
        1/12 --------------------------------------------------------------
                                          2
                                         h