Differentiation matrices - Fundamentals of Numerical Computation

In Finite differences we used finite differences to turn a discrete collection of function values into an estimate of the derivative of the function at a point. Just as with differentiation in elementary calculus, we can generalize differences at a point into an operation that maps discretized functions to discretized derivatives.

The derivative of a function is a unique result. The same is not true for our finite-dimensional approximation of the derivative, though.

10.3.1Matrices for finite differences¶

We first discretize the interval $x\in[a,b]$ into equal pieces of length $h=(b-a)/n$ , leading to the nodes

x_i=a+i h, \qquad i=0,\ldots,n.

(10.3.1)

Note again that our node indexing scheme starts at zero. Our goal is to find a vector $\mathbf{g}$ such that $g_i \approx f'(x_i)$ for $i=0,\ldots,n$ . Our first try is the forward difference formula (5.4.2),

g_i = \frac{f_{i+1}-f_i}{h}, \qquad t=0,\ldots,n-1.

(10.3.2)

However, this leaves $g_n$ undefined, because the formula would refer to the unavailable value $f_{n+1}$ . For $g_n$ we could resort to the backward difference

g_n = \frac{f_{n}-f_{n-1}}{h}.

(10.3.3)

We can summarize the entire set of formulas by defining

\mathbf{f} = \begin{bmatrix} f(x_0) \\[1mm] f(x_1) \\[1mm] \vdots \\[1mm] f(x_{n-1}) \\[1mm] f(x_n) \end{bmatrix}, \quad

(10.3.4)

and then the vector equation

\begin{bmatrix} f'(x_0) \\[1mm] f'(x_1) \\[1mm] \vdots \\[1mm] f'(x_{n-1}) \\[1mm] f'(x_n) \end{bmatrix} \approx \mathbf{D}_x \mathbf{f}, \qquad \mathbf{D}_x = \frac{1}{h} \begin{bmatrix} -1 & 1 & & & \\[1mm] & -1 & 1 & & \\[1mm] & & \ddots & \ddots & \\[1mm] & & & -1 & 1 \\[1mm] & & & -1 & 1 \end{bmatrix}.

(10.3.5)

The matrix $\mathbf{D}_x$ in (10.3.5) is a finite-difference differentiation matrix. Here as elsewhere, elements of $\mathbf{D}_x$ that are not shown are zero. Each row of $\mathbf{D}_x$ gives the weights of a finite-difference formula being used at one of the nodes.

The differentiation matrix $\mathbf{D}_x$ in (10.3.5) is not a unique choice. In fact, it’s about the least accurate choice possible, as explained in Convergence of finite differences. We are theoretically free to use whatever finite-difference formulas we want in each row, such as those in Table 5.4.1 and Table 5.4.2, although it makes sense to choose rows that are as similar as possible. For instance, using second-order centered differences where possible and second-order one-sided formulas at the boundary points leads to

\mathbf{D}_x = \frac{1}{h} \begin{bmatrix} -\frac{3}{2} & 2 & -\frac{1}{2} & & & \\[1mm] -\frac{1}{2} & 0 & \frac{1}{2} & & & \\[1mm] & -\frac{1}{2} & 0 & \frac{1}{2} & & \\ & & \ddots & \ddots & \ddots & \\ & & & -\frac{1}{2} & 0 & \frac{1}{2} \\[1mm] & & & \frac{1}{2} & -2 & \frac{3}{2} \end{bmatrix}.

(10.3.6)

The differentiation matrices so far are all sparse and banded, i.e., all the nonzero values are along diagonals close to the main diagonal.

10.3.2Second derivative¶

Similarly, we can define differentiation matrices for second derivatives. For example,

\begin{bmatrix} f''(x_0) \\[1mm] f''(x_1) \\[1mm] f''(x_2) \\[1mm] \vdots \\[1mm] f''(x_{n-1}) \\[1mm] f''(x_n) \end{bmatrix} \approx \frac{1}{h^2} \begin{bmatrix} 2 & -5 & 4 & -1 & & \\[1mm] 1 & -2 & 1 & & & \\[1mm] & 1 & -2 & 1 & & \\[1mm] & & \ddots & \ddots & \ddots & \\[1mm] & & & 1 & -2 & 1 \\[1mm] & & -1 & 4 & -5 & 2 \end{bmatrix} \begin{bmatrix} f(x_0) \\[1mm] f(x_1) \\[1mm] f(x_2) \\[1mm] \vdots \\[1mm] f(x_{n-1}) \\[1mm] f(x_n) \end{bmatrix} = \mathbf{D}_{xx} \mathbf{f}.

(10.3.7)

We have multiple choices again for $\mathbf{D}_{xx}$ . Moreover, while it is possible to define $\mathbf{D}_{xx}$ as $\mathbf{D}_x^2$ for some first-derivative $\mathbf{D}_x$ , it is not required---and often, not desirable---to do so, because it may put nonzeros farther from the diagonal than is necessary (see Finite differences).

10.3.3Implementation¶

Together the matrices (10.3.6) and (10.3.7) give second-order approximations of the first and second derivatives at all nodes. These matrices, as well as the nodes $x_0,\ldots,x_n$ , are returned by Function 10.3.1.

Algorithm 10.3.1 (diffmat2)

Julia

MATLAB

Python

Second-order differentiation matrices

diffmat2.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
"""
    diffmat2(n, xspan)

Compute 2nd-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 diffmat2(n, xspan)
    a, b = xspan
    h = (b - a) / n
    x = [a + i * h for i in 0:n]   # nodes

    # Define most of Dₓ by its diagonals.
    dp = fill(0.5 / h, n)        # superdiagonal
    dm = fill(-0.5 / h, n)       # subdiagonal
    Dₓ = diagm(-1 => dm, 1 => dp)

    # Fix first and last rows.
    Dₓ[1, 1:3] = [-1.5, 2, -0.5] / h
    Dₓ[n+1, n-1:n+1] = [0.5, -2, 1.5] / h

    # Define most of Dₓₓ by its diagonals.
    d0 = fill(-2 / h^2, n + 1)    # main diagonal
    dp = ones(n) / h^2         # super- and subdiagonal
    Dₓₓ = diagm(-1 => dp, 0 => d0, 1 => dp)

    # Fix first and last rows.
    Dₓₓ[1, 1:4] = [2, -5, 4, -1] / h^2
    Dₓₓ[n+1, n-2:n+1] = [-1, 4, -5, 2] / h^2

    return x, Dₓ, Dₓₓ
end

Second-order differentiation matrices

diffmat2.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function [x,Dx,Dxx] = diffmat2(n,xspan)
%DIFFMAT2   Second-order accurate differentiation matrices.
% Input:
%   n       number of subintervals (one less than the number of nodes)
%   xspan   interval endpoints
% Output:
%   x       equispaced nodes
%   Dx      matrix for first derivative 
%   Dxx     matrix for second derivative 

a = xspan(1);  b = xspan(2);
h = (b-a)/n;
x = a + h*(0:n)';   % nodes

% Define most of Dx by its diagonals.
dp = 0.5*ones(n,1)/h;       % superdiagonal
dm = -0.5*ones(n,1)/h;      % subdiagonal
Dx = diag(dm,-1) + diag(dp,1);

% Fix first and last rows.
Dx(1,1:3) = [-1.5,2,-0.5]/h;
Dx(n+1,n-1:n+1) = [0.5,-2,1.5]/h;

% Define most of Dxx by its diagonals.
d0 =  -2*ones(n+1,1)/h^2;   % main diagonal
dp =  ones(n,1)/h^2;        % superdiagonal and subdiagonal
Dxx = diag(dp,-1) + diag(d0) + diag(dp,1);

% Fix first and last rows.
Dxx(1,1:4) = [2,-5,4,-1]/h^2;
Dxx(n+1,n-2:n+1) = [-1,4,-5,2]/h^2;

Second-order differentiation matrices

diffmat2.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def diffmat2(n, xspan):
    """
    diffmat2(n, xspan)

    Compute 2nd-order-accurate differentiation matrices on n+1 points in the
    interval xspan. Return a vector of nodes, and the matrices for the first
    and second derivatives.
    """
    a, b = xspan
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)  # nodes

    # Define most of Dx by its diagonals.
    dp = 0.5 / h * np.ones(n)  # superdiagonal
    dm = -0.5 / h * np.ones(n)  # subdiagonal
    Dx = np.diag(dm, -1) + np.diag(dp, 1)

    # Fix first and last rows.
    Dx[0, :3] = np.array([-1.5, 2, -0.5]) / h
    Dx[-1, -3:] = np.array([0.5, -2, 1.5]) / h

    # Define most of Dxx by its diagonals.
    d0 = -2 / h**2 * np.ones(n + 1)  # main diagonal
    dp = np.ones(n) / h**2  # superdiagonal and subdiagonal
    Dxx = np.diag(d0, 0) + np.diag(dp, -1) + np.diag(dp, 1)

    # Fix first and last rows.
    Dxx[0, :4] = np.array([2, -5, 4, -1]) / h**2
    Dxx[-1, -4:] = np.array([-1, 4, -5, 2]) / h**2

    return x, Dx, Dxx

Example 10.3.1 (Differentiation matrices)

Julia

MATLAB

Python

Example 10.3.1

We test first-order and second-order differentiation matrices for the function $x + \exp(\sin 4x)$ over $[-1,1]$ .

f = x -> x + exp(sin(4 * x));

For reference, here are the exact first and second derivatives.

df_dx = x -> 1 + 4 * exp(sin(4x)) * cos(4x);
d2f_dx2 = x -> 4 * exp(sin(4x)) * (4 * cos(4x)^2 - 4 * sin(4x));

We discretize on equally spaced nodes and evaluate $f$ at the nodes.

t, Dₓ, Dₓₓ = FNC.diffmat2(18, [-1, 1])
y = f.(t);

Then the first two derivatives of $f$ each require one matrix-vector multiplication.

yₓ = Dₓ * y
yₓₓ = Dₓₓ * y;

The results show poor accuracy for this small value of $n$ .

using Plots
plot(df_dx, -1, 1, layout = 2, xaxis = (L"x"), yaxis = (L"f'(x)"))
scatter!(t, yₓ, subplot = 1)
plot!(d2f_dx2, -1, 1, subplot = 2, xaxis = (L"x"), yaxis = (L"f''(x)"))
scatter!(t, yₓₓ, subplot = 2)

A convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.

n = @. round(Int, 2^(4:0.5:11))
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
    t, Dₓ, Dₓₓ = FNC.diffmat2(n, [-1, 1])
    y = f.(t)
    err[k, 1] = norm(df_dx.(t) - Dₓ * y, Inf)
    err[k, 2] = norm(d2f_dx2.(t) - Dₓₓ * y, Inf)
end
plot(n, err, m = :o, label = [L"f'" L"f''"])
plot!(n, 10 * 10 * n .^ (-2);
    l = (:dash, :black),
    label = "2nd order",
    xaxis = (:log10, "n"),
    yaxis = (:log10, "max error"),
    title = "Convergence of finite differences")

Example 10.3.1

We test first-order and second-order differentiation matrices for the function $x + \exp(\sin 4x)$ over $[-1,1]$ .

f = @(x) x + exp(sin(4*x));

For reference, here are the exact first and second derivatives.

df_dx = @(x) 1 + 4 * exp(sin(4*x)) .* cos(4*x);
d2f_dx2 = @(x) 4 * exp(sin(4*x)) .* (4*cos(4*x).^2 - 4*sin(4*x));

We discretize on equally spaced nodes and evaluate $f$ at the nodes.

[t, Dx, Dxx] = diffmat2(12, [-1 1]);
y = f(t);

Then the first two derivatives of $f$ each require one matrix-vector multiplication.

yx = Dx * y;
yxx = Dxx * y;

The results show poor accuracy for this small value of $n$ .

clf,  subplot(2, 1, 1)
fplot(df_dx, [-1, 1]),  hold on
plot(t, yx, 'ko')
xlabel('x'),  ylabel('f''(x)')
subplot(2, 1, 2)
fplot(d2f_dx2, [-1, 1]),  hold on
plot(t, yxx, 'ko')
xlabel('x'),  ylabel('f''''(x)')

An convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.

n = round( 2 .^ (4:.5:11)' );
err = zeros(length(n), 2);
for k = 1:length(n)
    [t, Dx, Dxx] = diffmat2(n(k), [-1, 1]);
    y = f(t);
    err(k, 1) = norm(df_dx(t) - Dx * y, Inf);
    err(k, 2) = norm(d2f_dx2(t) - Dxx * y, Inf);
end
clf
loglog(n, err, 'o-'), hold on
loglog(n, 100 * n.^(-2), 'k--')
legend("f'", "f''", '2nd order')
xlabel('n'),  ylabel('max error')
title('Convergence of finite differences')

Example 10.3.1

We test first-order and second-order differentiation matrices for the function $x + \exp(\sin 4x)$ over $[-1,1]$ .

f = lambda x: x + exp( sin(4 * x) )

For reference, here are the exact first and second derivatives.

df_dx = lambda x: 1 + 4 * exp(sin(4 * x)) * cos(4 * x)
d2f_dx2 = lambda x: 4 * exp(sin(4 * x)) * (4 * cos(4 * x)**2 - 4 * sin(4 * x))

We discretize on equally spaced nodes and evaluate $f$ at the nodes.

t, Dx, Dxx = FNC.diffmat2(12, [-1, 1])
y = f(t)

Then the first two derivatives of $f$ each require one matrix-vector multiplication.

yx = Dx @ y
yxx = Dxx @ y

The results show poor accuracy for this small value of $n$ .

x = linspace(-1, 1, 500)
subplot(2, 1, 1)
plot(x, df_dx(x))
plot(t, yx, "ko")
xlabel("$x$"),  ylabel("$f'(x)$")

subplot(2, 1, 2)
plot(x, d2f_dx2(x))
plot(t, yxx, "ko")
xlabel("$x$"),  ylabel("$f''(x)$");

A convergence experiment confirms the order of accuracy. Because we expect an algebraic convergence rate, we use a log-log plot of the errors.

N = array([int(2**k) for k in arange(4, 11.5, 0.5)])
err1 = zeros(len(N))
err2 = zeros(len(N))
for k, n in enumerate(N):
    t, Dx, Dxx = FNC.diffmat2(n, [-1, 1])
    y = f(t)
    err1[k] = norm(df_dx(t) - Dx @ y, inf)
    err2[k] = norm(d2f_dx2(t) - Dxx @ y, inf)

loglog(N, err1, "-o", label="$f'$")
loglog(N, err2, "-o", label="$f''$")
plot(N, 10 * 10 / N**2, "k--", label="2nd order")
xlabel("$n$"),  ylabel("max error")
legend(loc="lower left")
title("Convergence of finite differences");

10.3.4Spectral differentiation¶

Recall that finite-difference formulas are derived in three steps:

Choose a node index set $S$ near node $i$ .
Interpolate with a polynomial using the nodes in $S$ .
Differentiate the interpolant and evaluate at node $i$ .

We can modify this process by using a global interpolant, either polynomial or trigonometric, as in Chapter 9. Rather than choosing a different index set for each node, we use all of the nodes each time.

In a nonperiodic setting we use Chebyshev second-kind points for stability:

x_k = -\cos\left(\frac{k \pi}{n}\right), \qquad k=0,\ldots,n;

(10.3.8)

then the resulting Chebyshev differentiation matrix has entries

\begin{gathered} D_{00} = \dfrac{2n^2+1}{6}, \qquad D_{n n} = -\dfrac{2n^2+1}{6}, \\ D_{ij} = \begin{cases} -\dfrac{x_i}{2(1-x_i^2)}, & i=j, \\[4mm] \dfrac{c_i}{c_j}\, \dfrac{(-1)^{i+j}}{x_i-x_j}, & i\neq j, \end{cases} \end{gathered}

(10.3.9)

where $c_0=c_n=2$ and $c_i=1$ for $i=1,\ldots,n-1$ .

The simplest way to compute a second derivative is by via $\mathbf{D}_{xx} = \mathbf{D}_x^2$ , as there is no longer any concern about the sparsity of the result.

Function 10.3.2 returns these two matrices. The function uses a change of variable to transplant the standard $[-1,1]$ for Chebyshev nodes to any $[a,b]$ . It also takes a different approach to computing the diagonal elements of $\mathbf{D}_x$ than the formulas in (10.3.9) (see Exercise 5).

Algorithm 10.3.2 (diffcheb)

Julia

MATLAB

Python

Chebyshev differentiation matrices

diffcheb.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
"""
    diffcheb(n, xspan)

Compute Chebyshev 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 diffcheb(n, xspan)
    x = [-cos(k * π / n) for k in 0:n]    # nodes in [-1,1]

    # Off-diagonal entries.
    c = [2; ones(n - 1); 2]    # endpoint factors
    dij = (i, j) -> (-1)^(i + j) * c[i+1] / (c[j+1] * (x[i+1] - x[j+1]))
    Dₓ = [dij(i, j) for i in 0:n, j in 0:n]

    # Diagonal entries.
    Dₓ[isinf.(Dₓ)] .= 0         # fix divisions by zero on diagonal
    s = sum(Dₓ, dims = 2)
    Dₓ -= diagm(s[:, 1])         # "negative sum trick"

    # Transplant to [a,b].
    a, b = xspan
    x = @. a + (b - a) * (x + 1) / 2
    Dₓ = 2 * Dₓ / (b - a)             # chain rule

    # Second derivative.
    Dₓₓ = Dₓ^2
    return x, Dₓ, Dₓₓ
end

Chebyshev differentiation matrices

diffcheb.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
function [x,Dx,Dxx] = diffcheb(n,xspan)
%DIFFCHEB   Chebyshev differentiation matrices.
% Input:
%   n      number of subintervals (integer)
%   xspan  interval endpoints (vector)
% Output:
%   x      Chebyshev nodes in domain (length n+1)
%   Dx     matrix for first derivative (n+1 by n+1)
%   Dxx    matrix for second derivative (n+1 by n+1)

x = -cos( (0:n)'*pi/n );    % nodes in [-1,1]
Dx = zeros(n+1);
c = [2; ones(n-1,1); 2];    % endpoint factors
i = (0:n)';                 % row indices

% Off-diagonal entries
for j = 0:n
  num = c(i+1).*(-1).^(i+j);
  den = c(j+1)*(x-x(j+1));
  Dx(:,j+1) = num./den;
end

% Diagonal entries
Dx(isinf(Dx)) = 0;          % fix divisions by zero on diagonal
Dx = Dx - diag(sum(Dx,2));  % "negative sum trick" 
  
% Transplant to [a,b]
a = xspan(1);  b = xspan(2);
x = a + (b-a)*(x+1)/2;
Dx = 2*Dx/(b-a);

% Second derivative
Dxx = Dx^2;

Chebyshev differentiation matrices

diffcheb.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def diffcheb(n, xspan):
    """
    diffcheb(n, xspan)

    Compute Chebyshev differentiation matrices on n+1 points in the
    interval xspan. Return a vector of nodes, and the matrices for the first
    and second derivatives.
    """
    x = -np.cos(np.arange(n + 1) * np.pi / n)  # nodes in [-1,1]
    Dx = np.zeros([n + 1, n + 1])
    c = np.hstack([2.0, np.ones(n - 1), 2.0])  # endpoint factors

    # Off-diagonal entries
    Dx = np.zeros([n + 1, n + 1])
    for i in range(n + 1):
        for j in range(n + 1):
            if i != j:
                Dx[i, j] = (-1) ** (i + j) * c[i] / (c[j] * (x[i] - x[j]))

    # Diagonal entries by the "negative sum trick"
    for i in range(n + 1):
        Dx[i, i] = -np.sum([Dx[i, j] for j in range(n + 1) if j != i])

    # Transplant to [a,b]
    a, b = xspan
    x = a + (b - a) * (x + 1) / 2
    Dx = 2 * Dx / (b - a)

    # Second derivative
    Dxx = Dx @ Dx

    return x, Dx, Dxx

Example 10.3.2 (Chebyshev differentiation matrices)

Julia

MATLAB

Python

Example 10.3.2

Here is a $4\times 4$ Chebyshev differentiation matrix.

t, Dₓ = FNC.diffcheb(3, [-1, 1])
Dₓ

4×4 Matrix{Float64}:
 -3.16667    4.0       -1.33333    0.5
 -1.0        0.333333   1.0       -0.333333
  0.333333  -1.0       -0.333333   1.0
 -0.5        1.33333   -4.0        3.16667

We again test the convergence rate.

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

n = 5:5:70
err1 = zeros(size(n))
err2 = zeros(size(n))
for (k, n) in enumerate(n)
    t, Dₓ, Dₓₓ = FNC.diffcheb(n, [-1, 1])
    y = f.(t)
    err1[k] = norm(df_dx.(t) - Dₓ * y, Inf)
    err2[k] = norm(d2f_dx2.(t) - Dₓₓ * y, Inf)
end

Since we expect a spectral convergence rate, we use a semi-log plot for the error.

plot(n, [err1 err2]; m = :o,
    label=[L"f'" L"f''"],
    xaxis=(L"n"),  yaxis = (:log10, "max error"),
    title="Convergence of Chebyshev derivatives")

Example 10.3.2

Here is a $4\times 4$ Chebyshev differentiation matrix.

[t, Dx] = diffcheb(3, [-1, 1]);
format rat
Dx

We again test the convergence rate.

f = @(x) x + exp(sin(4*x));
df_dx = @(x) 1 + 4 * exp(sin(4*x)) .* cos(4*x);
d2f_dx2 = @(x) 4 * exp(sin(4*x)) .* (4*cos(4*x).^2 - 4*sin(4*x));

n = 5:5:70;
err = zeros(length(n), 2);
for k = 1:length(n)
    [t, Dx, Dxx] = diffcheb(n(k), [-1, 1]);
    y = f(t);
    err(k, 1) = norm(df_dx(t) - Dx * y, Inf);
    err(k, 2) = norm(d2f_dx2(t) - Dxx * y, Inf);
end

Since we expect a spectral convergence rate, we use a semi-log plot for the error.

clf,  format
semilogy(n, err, 'o-'), hold on
legend("f'", "f''")
xlabel('n'),  ylabel('max error')
title('Convergence of finite differences')

Example 10.3.2

Here is a $4\times 4$ Chebyshev differentiation matrix.

t, Dx, Dxx = FNC.diffcheb(3, [-1, 1])
print(Dx)

[[-3.16666667  4.         -1.33333333  0.5       ]
 [-1.          0.33333333  1.         -0.33333333]
 [ 0.33333333 -1.         -0.33333333  1.        ]
 [-0.5         1.33333333 -4.          3.16666667]]

We again test the convergence rate.

f = lambda x: x + exp(sin(4 * x))
df_dx = lambda x: 1 + 4 * exp(sin(4 * x)) * cos(4 * x)
d2f_dx2 = lambda x: 4 * exp(sin(4 * x)) * (4 * cos(4 * x) ** 2 - 4 * sin(4 * x))

N = range(5, 75, 5)
err1 = zeros(len(N))
err2 = zeros(len(N))
err = zeros((len(N), 2))
for k, n in enumerate(N):
    t, Dx, Dxx = FNC.diffcheb(n, [-1, 1])
    y = f(t)
    err[k, 0] = norm(df_dx(t) - Dx @ y, inf)
    err[k, 1] = norm(d2f_dx2(t) - Dxx @ y, inf)

Since we expect a spectral convergence rate, we use a semi-log plot for the error.

semilogy(N, err, "-o")
xlabel("$n$"), ylabel("max error")
legend(["$f'$", "$f''$"], loc="lower left")
title("Convergence of Chebyshev derivatives");

According to Theorem 9.3.1, the convergence of polynomial interpolation to $f$ using Chebyshev nodes is spectral if $f$ is analytic (at least having infinitely many derivatives) on the interval. The derivatives of $f$ are also approximated with spectral accuracy.

10.3.5Exercises¶

(a) ✍ Calculate $\mathbf{D}_x^2$ using (10.3.5) to define $\mathbf{D}_x$ .
(b) ⌨ Repeat the experiment of Demo 10.3.1, but using this version of $\mathbf{D}_x^2$ to estimate $f''$ . What is the apparent order of accuracy in $f''$ ?
(a) ✍ Find the derivative of $f(x) =\operatorname{sign}(x)x^2$ on the interval $[-1,1]$ . (If this gives you trouble, use an equivalent piecewise definition of $f$ .) What is special about this function at $x=0$ ?
(b) ⌨ Adapt Demo 10.3.1 to operate on the function from part (a), computing only the first derivative. What is the observed order of accuracy?
(c) ✍ Show that for even values of $n$ , there is only one node at which the error for computing $f'$ in part (b) is nonzero.
⌨ 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.
(a) Modify Function 10.3.1 to a function diffmat4, which outputs fourth-order accurate differentiation matrices. You may want to use Function 5.4.1.
(b) Repeat the experiment of Demo 10.3.1 using diffmat4 in place of Function 10.3.1, and compare observed errors to fourth-order accuracy.
✍ Explain in detail how lines 23-24 in Function 10.3.2 correctly change the interval from $[-1,1]$ to $[a,b]$ .

(a) ✍ What is the derivative of a constant function?
(b) ✍ Explain why for any reasonable differentiation matrix $\mathbf{D}$ , we should find $\displaystyle \sum_{j=0}^nD_{ij}=0$ for all $i$ .
(c) ✍ What does this have to do with Function 10.3.2? Refer to specific line(s) in the function for your answer.
Define the $(n+1)\times (n+1)$ matrix $\mathbf{T} = \displaystyle \begin{bmatrix} 1 & & & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & 1 & \cdots & 1 \end{bmatrix}$ .
(a) ✍ Write out $\mathbf{T}\mathbf{u}$ for a generic vector $\mathbf{u}$ (start with a zero index). How is this like integration?
(b) ✍ Find the inverse of $\mathbf{T}$ for any $n$ . (You can use Julia to find the pattern, but show that the result is correct in general.) What does this have to do with the inverse of integration?

Preface

Shooting

Preface

Collocation for linear problems