Error in result of finite-difference approximation when refining

Question

I have calculated the first derivative of following equation using Euler method (first order), Three point Finite Difference method (second order) and Four point Finite Difference method (third order). $$f(x)=e^{-32*(x-5)^2}$$ in the domain $0 \leq x \leq 10$.

I discretized the domain by 50, 500, 5000 and 50000 number of points and calculated the solution. Then I compared the solution with analytical solution and calculated the maximum absolute difference in analytical and numerical result in that region and plotted the graph shown below (sorry for bad quality of graph). here y-axis is log(error) and x-axis is log(number of grid points).

My question is:

Higher order methods are supposed to be more accurate than lower order methods. In the graph error in third order method is more than second and first order method for course grid (50 and 500 points) and third order method is more accurate than second and first order method in fine grid ( 5000 and 50000) points). Please explain me why this is happening?

thank you!

Answer 1

High order methods are not necessarily always more accurate than low order methods; they simply have a more rapid convergence rate. For example, if you take the Taylor expansion of a function at a point $x$ and evaluate it at $x+h$, you have

$$f(x+h) = f(x) + f'(x)h + \frac{f''(x)}{2}h^2 + \ldots$$

The order of a finite difference method is determined by the power of $h$ that is eliminated using finite difference coefficients; for example, rearranging gives

$$\left|\frac{f(x+h)-f(x)}{h} - f'(x)\right| = \left| \frac{f''(x)}{2}h + \frac{f'''(x)}{6}h^2 + \ldots\right|$$

The left term is error ,and you can make the asymptotic argument that if $h \ll 1$, then $h \gg h^2$, so that for $h$ small enough,

$$\left|\frac{f(x+h)-f(x)}{h} - f'(x)\right| \approx \left| \frac{f''(x)}{2}\right| h.$$

This truncation error term implies that error decreases linearly with grid size ($h$), or that the error is $O(h)$. High order methods give that the error is $O(h^r)$ for some order $r$; again, this doesn't imply that high order method errors are smaller, just that the error decreases faster with $h$ than for a low order method.The magnitude of the truncation error depends specifically on magnitude of $ f'', f'''$ etc. we truncated, and power of $h$..

Problem you have considered here has complex behavior in derivatives. If your grid size is not small enough to nullify the effect of magnitude of derivative term you truncated, will give high error. Because of high rate of convergence of higher order methods, your 4th order method error decreases with a faster slope than the lower order methods.

Answer 2

You've chosen a very specific example where the function in question has a single extremely narrow peak and is almost zero everywhere else in the domain. This is an example of a stiff problem where you must take a very tight discretization to resolve the peak. Finite difference approximations are based on polynomial approximations to a curve. In the case where the curve is not locally well-approximated by a polynomial (such as the several grid points that straddle the peak in the Gaussian function here) this can lead to oscillations in the approximating polynomial (see this answer) resulting in a poor estimate of the derivative. Refining the grid eventually leads to a situation where the function is well approximated by a polynomial on the scale of the finite difference stencil and accuracy improves.

As discussed in Jesse Chan's answer, low order methods may actually have lower magnitude errors in some cases. However, I would be skeptical of the accuracy of either method in that case and suggest that you simply need higher resolution.

I tried to replicate your results, but I found that higher order finite differences were actually more accurate (based on RMS error) even in this example. Were you using one sided differences? Here's my matlab script if you want to mess with it.

clear

y = @(x) exp(-32*(x-5).^2);
dydx = @(x) -(64*(x-5)).*y(x);

x=linspace(0,10,50);

%first order (one sided) difference
x1 = x(1:end-1);
dy1 = (y(x(2:end))-y(x(1:end-1)))/diff(x(1:2));

%second order centered difference
x2 = x(2:end-1);
dy2 = (y(x(3:end))-y(x(1:end-2)))/(2*diff(x(1:2)));

%fourth order centered difference
x3 = x(3:end-2);
dy3 = (-y(x(5:end))+8*y(x(4:end-1))-8*y(x(2:end-3))+y(x(1:end-4)))/(12*diff(x(1:2)));

xRef = linspace(0,10,1e5);
dyRef = dydx(xRef);

figure(1)
clf
hold on
plot( xRef, dyRef, 'k' )
plot( x1, dy1,'o')
plot( x2, dy2,'x' )
plot( x3, dy3,'^' )
legend('Analytical','First order','Second order','Fourth order')
title('Computed derivative values')

figure(2)
clf
hold on
plot( x1, abs(dy1-dydx(x1)),'o')
plot( x2, abs(dy2-dydx(x2)),'x' )
plot( x2, abs(dy2-dydx(x2)),'^' )
legend('First order','Second order','Fourth order')
title('Error in the derivative computation')

% rms errors:
sqrt( sum((dy1-dydx(x1)).^2 ) / numel(x1) ) %first order
sqrt( sum((dy2-dydx(x2)).^2 ) / numel(x2) ) %second order
sqrt( sum((dy3-dydx(x3)).^2 ) / numel(x3) ) %fourth order

Please note that this is a slightly different situation than that in the related question. That question was about ODE solvers (Euler and RK). The reason I edited your answer there was to highlight the fact that this case is not typical and that in almost all cases where RK4 is stable it will outperform Euler for a given step size. The ODE analogue to your example here is a stiff ODE where you are not in the stability region of explicit solvers. As I said there, if your ODE solver is unstable all bets are off as far as accuracy (just as all the finite difference methods give poor results on a coarse grid here). Of course a stable low-order implicit method will typically have less error than an unstable explicit method.

Even though you may have "less error" for course grids with the lower order finite differences, the actual errors in your example are huge. In that case, neither solution should be used and you have no choice but to refine the grid or use a different method entirely.