Numerical Differentiation Formulas

In Physics, to simulate physical system, we usually encounter ordinary or partial differential equations. Let us consider a function or a field f(x) which satisfies a specified differential equation. To solve the differential equation, we first approximate the function f(x). The approximation can be in a form of Fourier series, a polynomial and a discrete functions. Today I will only consider the discrete function. The function f(x) is approximated in an interval (a,b) by (N+1) discrete values given by {f(i)} where f(x) = f(i*h) , h = (b-a)/N, i=0,1,2 ... N. Then we need to approximate derivatives in the differential equation.

The numerical differentiations that I often used are central finite difference schemes:

1) d f(x)/dx = [f(x+h)-f(x-h)]/2h + O(h^2)
= [f(i+1)-f(i-1)]/2h (second order accurate)

2) d^2 f(x)/dx^2 = [f(x+h)-2f(x) + f(x-h)]/h^2 + O(h^2)
= [f(i+1)-2f(i)+f(i-1)]/h^2 (second order accurate)

d^2 f(x)/dx^2 = [-f(x+2h)+16f(x+h)-30f(x) + 16f(x-h)-f(x-2h)]/(12h^2) + O(h^4)
= [-f(i+2)+16f(i+1)-30f(i)+16f(i-1)-f(i-2)]/(12h^2)
(fourth order accurate)

3) d^3 f(x)/dx^3 = [f(x+2h)-2f(x+h) +2 f(x-h)-f(x-2h)]/(2h^3) + O(h^2)
= [f(i+2)-2f(i+1)+2f(i-1)-f(i-2)]/(2h^3) (second order accurate)

As an example, Let us compute derivatives for trigonometric function sin(Pi x) in the interval (0,1). In C programming language, can be written as

--- num_diff.c ---

/* Numerical Differentiation
* by I Wayan Sudiarta Oct 2007
*/

#include<math.h>
#include<stdio.h>

#define N 20
#define PI 3.141592653589793

int main(void)
{
int i;
double h, hsq, hcb;
double f[N+1], f1[N+1], f2[N+1], f2b[N+1], f3[N+1];

h = 1.0/((double)(N));

hsq = h*h;
hcb = h*hsq;

/* Compute f[i] */
for(i=0;i<=N;i++){
f[i] = sin(i*PI*h);
}

/* Compute derivatives */
for(i=2;i<(N-1);i++){
f1[i] = (f[i+1]-f[i-1])/(2*h);
f2[i] = (f[i+1]-2*f[i]+f[i-1])/(hsq);
f2b[i]= (-f[i+2]+16*f[i+1]-30*f[i]+16*f[i-1]-f[i-2])/(12*hsq);
f3[i] = (f[i+2]-2*f[i+1]+2*f[i-1]-f[i-2])/(2*hcb);
}


/* Print results */
printf("x f1(x) f2(x) f2b(x) f3(x) \n");
for(i=2;i<(N-1);i++){
printf("%e %e %e %e %e \n", i*h, f1[i], f2[i], f2b[i], f3[i]);
}
return 0;
}

Ref: CRC Standard Mathematical Tables, 27th Edition, 1974