Centered approximations to $v''(x)$

AMath 586, Spring Quarter 2019 at the University of Washington. For other notebooks, see Index.html or the Index of all notebooks on Github.

This illustrates the centered difference approximation to the second derivative $v''(x)$ for a smooth function of $x$. This approximation is used for the spatial discretization of $u_{xx}(x,t)$ in solving PDEs.

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

Suppose we are given a function $v(x)$ that is $2\pi$-periodic, i.e. we know the values for $0 \leq x \leq 2\pi$, and we wish to approximate the derivative $v\prime(x)$ at points in this interval.

To be specific, choose $N$ equally spaced points $x_j = jh$ with equal spacing $h = 2\pi / N$, for $j=1,~2,~\ldots,~N$. Then we know the discrete function values $v_j = v(x_j)$ and want to compute approximations $w_j\approx v\prime(x_j)$ or $z_j \approx v\prime\prime(x_j)$.

Finite differences

The centered difference approximation to $v''(x) = \frac 1 {h^2} (v(x-h) - 2v(x) + v(x+h)) + O(h^2)$ can be written as $z = A v$ (so that $z_j \approx v\prime\prime(x_j)$) with $$ A = \frac{1}{h^2} \left(\begin{array}[rrrrrr] -2&1\\ 1&-2&1\\ &1&-2&1\\ &&&\ddots\\ &&&1&-2&1\\ &&&&1&-2 \end{array}\right) $$

Python code

We can use the spdiags function from scipy.sparse to set up a sparse matrix.

In [2]:
from scipy.sparse import spdiags

m = 9
h = 1. / (m+1)
x = linspace(h,m*h,m)  # interior grid points
e = ones((m,1))
data = hstack([e, -2*e, e]).T
diags = [-1, 0, 1]
T = spdiags(data, diags, m, m)
A = (1./h**2) * T
In [3]:
print("With m = %i, h = %6.5f" % (m,h))
print("Interior grid points: \n", x)
print("\nA = \n", h**2 * A.todense())

With m = 9, h = 0.10000
Interior grid points: 
 [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]

A = 
 [[-2.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 1. -2.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  1. -2.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  1. -2.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  1. -2.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  1. -2.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1. -2.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  1. -2.  1.]
 [ 0.  0.  0.  0.  0.  0.  0.  1. -2.]]

Test this out

Choose a function $v(x)$ with $v(0)=v(1)=0$ and work out the derivatives:

In [4]:
v_0 = lambda x: sin(x*(x-1))
v_1 = lambda x: (2*x-1)*cos(x*(x-1))
v_2 = lambda x: -(2*x-1)**2 * sin(x*(x-1)) + 2*cos(x*(x-1))

# fine grid for plotting function and true derivative:
x_fine = linspace(0,1.,1000)  
v_fine = v_0(x_fine)
z_fine = v_2(x_fine)
In [5]:
def test_vxx(m):
    h = 1. / (m+1)
    x = linspace(h,m*h,m)  # interior grid points
    e = ones((m,1))
    data = hstack([e, -2*e, e]).T
    diags = [-1, 0, 1]
    T = spdiags(data, diags, m, m)
    A = (1./h**2) * T
    v = v_0(x)      # evaluate v(x_j)
    z = A.dot(v)   # sparse matrix vector product A * v

    figure(figsize=(15,5))
    subplot(121)
    plot(x_fine,v_2(x_fine))
    plot(x,z,'ro')
    title("v''(x) and approximation z")

    subplot(122)
    error = z - v_2(x)
    plot(x,error,'ro-')
    title('h = %g, Max error = %g' % (h,norm(error,inf)))

test_vxx(19)
test_vxx(39)

Note that reducing $h$ by a factor of 2 reduces the error by roughly a factor of 4, as expected since the approximation should be second order accurate.

If you tried this on a function that doesn't satisfy $v(0)=v(1)=0$ there would be large errors at the the boundaries unless you adjusted $Av$ to incorporate the boundary values properly.

Periodic boundary conditions

If we assume periodic boundary conditions, the matrix is $$ A = \frac{1}{h^2} \left(\begin{array}[rrrrrr] -2&1&&&&1\\ 1&-2&1\\ &1&-2&1\\ &&&\ddots\\ &&&1&-2&1\\ 1&&&&1&-2 \end{array}\right) $$

Now the grid has to include one boundary point or the other (not both). Below we include the right boundary point.

In [6]:
m = 9
h = 1./(m+1)
x = linspace(h,1,m+1)  # interior grid points and right boundary

e = ones((m+1,1))
data = hstack([e, e, -2*e, e, e]).T
diags = [-m, -1, 0, 1, m]
T = spdiags(data, diags, m+1, m+1)
A_periodic = (1./h**2) * T

For small N we can print the D2 matrix without the $1/h^2$ factor so we can verify we did it right:

In [7]:
print("With m = %i, h = %6.5f" % (m,h))
print("There are %i grid points: \n" % (m+1), x)
print("\nA_periodic = \n", h**2 * A_periodic.todense())

With m = 9, h = 0.10000
There are 10 grid points: 
 [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]

A_periodic = 
 [[-2.  1.  0.  0.  0.  0.  0.  0.  0.  1.]
 [ 1. -2.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  1. -2.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1. -2.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  1. -2.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1. -2.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  1. -2.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1. -2.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  1. -2.  1.]
 [ 1.  0.  0.  0.  0.  0.  0.  0.  1. -2.]]

Test this out.

Choose a function $v(x)$ with a smooth periodic extension and compute exact derivatives:

In [8]:
v_0 = lambda x: exp(sin(2*pi*x))
v_1 = lambda x: 2*pi*cos(2*pi*x) * v_0(x)
v_2 = lambda x: ((2*pi*cos(2*pi*x))**2 - 4*pi**2*sin(2*pi*x)) * v_0(x)

# fine grid for plotting function and true derivative:
x_fine = linspace(0,1.,1000)  
v_fine = v_0(x_fine)
z_fine = v_2(x_fine)
In [9]:
def test_vxx(m):
    h = 1./(m+1)
    x = linspace(h,1,m+1)  # interior grid points and right boundary

    e = ones((m+1,1))
    data = hstack([e, e, -2*e, e, e]).T
    diags = [-m, -1, 0, 1, m]
    T = spdiags(data, diags, m+1, m+1)
    A_periodic = (1./h**2) * T

    v = v_0(x)      # evaluate v(x_j)
    z = A_periodic.dot(v)   # sparse matrix vector product A * v

    figure(figsize=(15,5))
    subplot(121)
    plot(x_fine,v_2(x_fine))
    plot(x,z,'ro')
    title("v''(x) and approximation z")

    subplot(122)
    error = z - v_2(x)
    plot(x,error,'ro-')
    title('h = %g, Max error = %g' % (h,norm(error,inf)))

test_vxx(19)
test_vxx(39)
test_vxx(79)

Again note that the error goes down by roughly a factor of 4 each time we double $m$.

In [ ]: