p01: Convergence of fourth order finite differences

Contents

p01: Convergence of fourth order finite differences#

Let us compute the derivative of

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

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) \]

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) \]

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'| \]
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');
_images/e41baccb8a07c1d86f7253f2035317c6fe393e08357da4cf9362b5117f527ee5.svg

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.

Contents