sandbox/Antoonvh/pade.c
8th-order accurate finite-differencing using a 5-point stencil
Solving the system
\displaystyle \beta f'_{-2} + \alpha f'_{-1} + f'+\alpha f'_{1} + \beta f'_{2} = a\frac{f_{1} - f_{-1}}{2 \Delta }+ b\frac{f_{2} - f_{-2}}{4\Delta},
with coefficients,
\displaystyle \alpha = \frac{4}{9}, \displaystyle \beta = \frac{1}{36}, \displaystyle a = \frac{40}{27}, \displaystyle b = \frac{25}{54},
yields an 8-th order accurate estimation for f' (Lele, 1991).
set logscale xy
set grid
set size square
set xr [4:1024]
set xlabel 'N'
set ylabel 'error'
plot 'out' u 1:2 t 'L1',\
'' u 1:3 t 'Max' ,\
1e10*x**-8 t '8th order'
reset
set logscale x
set grid
set size square
set xr [4:1024]
set yr [0:9]
set xlabel 'N'
set ylabel 'Iterations'
set key top left
plot 'out' u 1:4 t 'Cycles', '' u 1:5 t 'sweeps'
#define BGHOSTS 2
#include "grid/multigrid1D.h"
#include "poisson.h"
#define DF1(f) ((f[1] - f[-1])/(2*Delta))
#define DF2(f) ((f[2] - f[-2])/(4*Delta))
double alphaP = 4./9., betaP = 1./36.;
double aP = 40./27., bP = 25./54.;
double residual_pade (scalar * al, scalar * bl,
scalar * resl, void * data) {
scalar a = al[0], b = bl[0], res = resl[0];
double mr = 0;
foreach(reduction (+:mr)) {
res[] = (a[] + alphaP*(a[1] + a[-1]) + betaP*(a[2] + a[-2]) -
(aP*DF1(b) + bP*DF2(b)));
if (fabs(res[]) > mr)
mr = fabs(res[]);
}
return mr;
}
void relax_pade (scalar * al, scalar * resl,
int depth, void * data) {
scalar a = al[0], res = resl[0];
foreach_level(depth) {
a[] = -(res[] + alphaP*(a[-1] + a[1]) + betaP*(a[-2] + a[2]));
}
}
Set-up a test case
Differentiating a Gaussian pulse.
#define F(x) (exp(-sq(x)))
#define DF(x) (-2.*x*exp(-sq(x)))
scalar f[], dfdx[];
int main() {
L0 = 20;
X0 = -L0/2.;
for (N = 8; N <= 512; N *= 2) {
init_grid(N);
foreach()
f[] = F(x);
boundary ({f});
foreach()
dfdx[] = 4./3.*DF1(f) - DF2(f)/3.; //Guess
boundary ({dfdx}); // BC for the derivative...
mgstats mgs = mg_solve ({dfdx}, {f},
residual_pade, relax_pade);
double e = 0, em = 0;
foreach() {
double el = fabs(DF(x) - dfdx[]);
e += Delta*el;
if (el > em)
em = el;
}
printf ("%d %g %g %d %d\n", N, e, em, mgs.i, mgs.nrelax);
TOLERANCE = 1e-3*em;
}
}
Note
A 6th order accurate approximation for f' can be obtained with the coefficients.
\displaystyle a = \frac{2}{3}(\alpha + 2), \displaystyle b = \frac{1}{3}(4\alpha - 1), \displaystyle \alpha = \frac{1}{3}, \displaystyle \beta = 0.
Which has a 3-point relaxation stencil.
Futher, a 4th order method is achievable on a 3-point stencil.
Reference
S.K. Lele, Compact Finite Difference Schemes with Spectral-like Resolution (1991), in JCP.