p01: Convergence of fourth order finite differences

Let us compute the derivative of

u(x) = \exp(\sin(x)), \qquad x \in [-\pi,\pi]

(1)

using fourth order finite difference scheme

u'(x_j) \approx w_j = \frac{1}{h} \left( \frac{1}{12} u_{j-2} - \frac{2}{3} u_{j-1} + \frac{2}{3} u_{j+1} - \frac{1}{12} u_{j+2} \right)

(2)

using periodic boundary conditions.

#%config InlineBackend.figure_format='svg'
from scipy.sparse import coo_matrix
from numpy import arange,pi,exp,sin,cos,ones,inf
from numpy.linalg import norm
from matplotlib.pyplot import figure,loglog,semilogy,text,grid,xlabel,ylabel,title

The differentiation matrix is anti-symmetric; we first construct the upper triangular part $U$ and get

D = \frac{1}{h} ( U - U^\top)

(3)

The matrix is stored in coo_format which is a sparse matrix format.

The error in the derivative approximation is measured as

E_N = \max_{1 \le j \le N} |w_j - u_j'|

(4)

Nvec = 2**arange(3,13)
for N in Nvec:
    h = 2*pi/N
    x = -pi + arange(1,N+1)*h
    u = exp(sin(x))   # function values
    uprime = cos(x)*u # Exact derivative
    e = ones(N)
    e1 = arange(0,N)
    e2 = arange(1,N+1); e2[N-1]=0
    e3 = arange(2,N+2); e3[N-2]=0; e3[N-1]=1;
    D = coo_matrix((2*e/3,(e1,e2)),shape=(N,N)) \
        - coo_matrix((e/12,(e1,e3)),shape=(N,N))
    D = (D - D.T)/h
    error = norm(D@u-uprime,inf)
    loglog(N,error,'or')
    
loglog(Nvec,100*Nvec**(-4.0),'--')
text(105,5e-8,'$N^{-4}$')
grid(True); xlabel('N'); ylabel('error')
title('Convergence of fourth-order finite difference');

With about 4000 points, the error is of order 10^-12.

Exercise¶

Find sixth and eigth order finite differences, implement them and make error plot as above.

Exercise¶

Construct the differentiation matrix using scipy.toeplitz function. You need to pass the first column and first row.