/** # 8th-order accurate finite-differencing using a 5-point stencil Solving the system $$\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, $$\alpha = \frac{4}{9},$$ $$\beta = \frac{1}{36},$$ $$a = \frac{40}{27},$$ $$b = \frac{25}{54},$$ yields an 8-th order accurate estimation for $f'$ (Lele, 1991). ~~~gnuplot Convergence is OK 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' ~~~ ~~~gnuplot Convergence is OK 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. $$a = \frac{2}{3}(\alpha + 2),$$ $$b = \frac{1}{3}(4\alpha - 1),$$ $$\alpha = \frac{1}{3},$$ $$\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. */