Finite Difference Stencil Error Tests

AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.

These notebooks are all available on Github.

This notebook demonstrates how the error behaves as $h = \Delta x$ is reduced in two different approximations to $u''(\bar x)$, the standard second-order accurate 3-point approximation and a fourth-order accurate centered 5-point approximation.

A table of errors is generated and then log-log plots of the error are compared with the behavior expected from the general expression for the truncation error of each method.

In [1]:
%matplotlib inline
In [2]:
from pylab import *

import fdstencil function

The python module fdstencil_functions.py contains the functions that are defined in the notebook fdstencil.html, and was generated by exporting that notebook as Python and then editing out everything except the function definitions:

In [3]:
from fdstencil_functions import fdstencil

Test function and its derivatives:

In [4]:
u = lambda x: sin(2*x)
u2 = lambda x: -4*sin(2*x)
u4 = lambda x: 16*sin(2*x)
u6 = lambda x: -64*sin(2*x)

Test centered 3-point approximation

In [5]:
x0 = 1.
u2_true = u2(x0)
h0 = 0.1
stencil_points = array([-1, 0, 1])
c = fdstencil(2, 0, stencil_points)
h_vals = []
errors = []
print('\n    h                 error ')
for n in range(15):
    h = h0 / 2**n
    h_vals.append(h)
    x = x0 + h*stencil_points
    u2_approx = dot(c,u(x)) / h**2
    error = abs(u2_approx - u2_true)
    errors.append(error)
    print('%10.8f   %20.16f' % (h,error))

Stencil for approximation to order 2 derivative at U_{0} is:
[ + 1 U_{-1} - 2 U_{0} + 1 U_{1}] / h^2

    h                 error 
0.10000000     0.0121078119449098
0.05000000     0.0030299812726771
0.02500000     0.0007576847129194
0.01250000     0.0001894330170877
0.00625000     0.0000473589895229
0.00312500     0.0000118397966027
0.00156250     0.0000029600019253
0.00078125     0.0000007399708331
0.00039063     0.0000001844514692
0.00019531     0.0000000476634661
0.00009766     0.0000000214700187
0.00004883     0.0000000367376423
0.00002441     0.0000002695682859
0.00001221     0.0000002892252589
0.00000610     0.0000010342833185
In [6]:
loglog(h_vals, errors, 'bx-', label='Observed errors')
xlim(1e-6,1)
ylim(1e-17,1)
grid(True)
xlabel('h = Delta x')
ylabel('abs(error)')

eref = 1/12. * array(h_vals)**2 * abs(u4(x0))
loglog(h_vals, eref, 'r-', label='Expected errors')
legend(loc='lower right')
Out[6]:
<matplotlib.legend.Legend at 0x819153dd8>

Test centered 5-point approximation

In [7]:
x0 = 1.
u2_true = u2(x0)
h0 = 0.1
stencil_points = array([-2, -1, 0, 1, 2])
c = fdstencil(2, 0, stencil_points)
h_vals = []
errors = []
print('\n    h                 error ')
for n in range(10):
    h = h0 / 2**n
    h_vals.append(h)
    x = x0 + h*stencil_points
    u2_approx = dot(c,u(x)) / h**2
    error = abs(u2_approx - u2_true)
    errors.append(error)
    print('%10.8f   %20.16f' % (h,error))

Stencil for approximation to order 2 derivative at U_{0} is:
[ - 0.0833333 U_{-2} + 1.33333 U_{-1} - 2.5 U_{0} + 1.33333 U_{1} - 0.0833333 U_{2}] / h^2

    h                 error 
0.10000000     0.0000644306482318
0.05000000     0.0000040377153994
0.02500000     0.0000002525267924
0.01250000     0.0000000157885189
0.00625000     0.0000000009893348
0.00312500     0.0000000000968932
0.00156250     0.0000000001878426
0.00078125     0.0000000007335395
0.00039063     0.0000000010973373
0.00019531     0.0000000098284865
In [8]:
loglog(h_vals, errors, 'bx-', label='Observed errors')
xlim(1e-5,1)
ylim(1e-17,1)
grid(True)
xlabel('h = Delta x')
ylabel('abs(error)')

eref = 2/135. * array(h_vals)**4 * abs(u6(x0))
loglog(h_vals, eref, 'r-', label='Expected errors')
legend(loc='lower right')
Out[8]:
<matplotlib.legend.Legend at 0x8192d1cf8>