Exercise 10.3.3

To get a fourth-order accurate version of $\mathbf{D}_x$, five points per row are needed, including two special rows at each boundary. For a fourth-order $\mathbf{D}_{xx}$, five symmetric points per row are needed for interior rows and six points are needed for the rows near a boundary.

1. modify diffmat2() to a function diffmat4, which outputs fourth-order accurate differentiation matrices. (may use fdweights())

"""
    diffmat4(n,xspan)

Compute 4th-order-accurate differentiation matrices on `n`+1 points
in the interval `xspan`. Returns a vector of nodes and the matrices
for the first and second derivatives.
"""
function diffmat4(n,xspan)
    a,b = xspan
    h = (b-a)/n
    x = [ a + i*h for i in 0:n ]   # nodes

    # Define most of Dx by its diagonals.
    int_pts = FNC.fdweights(x[1:5].-x[3], 1)
    Dx = diagm(-2=>fill(int_pts[1],n-1),
               -1=>fill(int_pts[2],n),
               0=>fill(int_pts[3],n+1),
               1=>fill(int_pts[4],n),
               2=>fill(int_pts[5],n-1))
    
    # Boundaries 1, 2, n, n+1
    for i in 1:2
      Dx[i,1:5] = FNC.fdweights(x[1:5].-x[i],1)
      Dx[(n+i-1),((n+1)-4):(n+1)] =
        FNC.fdweights(x[((n+1)-4):(n+1)].-x[(n+i-1)],1)
    end
    
    # Define most of Dxx by its diagonals.
    int_pts = FNC.fdweights(x[1:5].-x[3], 2)
    Dxx = diagm(-2=>fill(int_pts[1],n-1),
               -1=>fill(int_pts[2],n),
               0=>fill(int_pts[3],n+1),
               1=>fill(int_pts[4],n),
               2=>fill(int_pts[5],n-1))

    # Fix first and last rows.
    for i in 1:2
      Dxx[i,1:6] = FNC.fdweights(x[1:6].-x[i],2)
      Dxx[(n+i-1),((n+1)-5):(n+1)] =
        FNC.fdweights(x[((n+1)-5):(n+1)].-x[(n+i-1)],2)
    end

    return x,Dx,Dxx
end

f = x -> x + exp(sin(4*x));
dfdx = x -> 1 + 4*exp(sin(4*x))*cos(4*x);
d2fdx2 = x ->  4*exp(sin(4*x))*(4*cos(4*x)^2-4*sin(4*x));

t, Dx, Dxx = diffmat4(18, [-1,1]);
y = f.(t);
Dx[1:5, 1:10]
Dxx[1:5, 1:10]

yx = Dx*y;
yxx = Dxx*y;
plot(dfdx,-1,1,layout=2,xaxis=(L"x"),yaxis=(L"f'(x)"));
scatter!(t, yx,subplot=1);
plot!(d2fdx2,-1,1,subplot=2,xaxis=(L"x"),yaxis=(L"f''(x)"));
scatter!(t, yxx,subplot=2)

n = @. round(Int,2^(4:.5:11) );
err1 = zeros(size(n));
err2 = zeros(size(n));
for (k,n) in enumerate(n)
    t,Dx,Dxx = diffmat4(n,[-1,1])
    y = f.(t)
    err1[k] = norm( dfdx.(t) - Dx*y, Inf )
    err2[k] = norm( d2fdx2.(t) - Dxx*y, Inf )
end
plot(n,[err1 err2],m=:o,label=[L"f'" L"f''"]);
plot!(n,10*10*n.^(-4),l=(:dash,:gray),label="4th order",
    xaxis=(:log10,"n"), yaxis=(:log10,"max error"),
    title="Convergence of finite differences")