Differentiation matrices

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.

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

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),

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

We can summarize the entire set of formulas by defining

and then the vector equation

Here as elsewhere, elements of $\mathbf{D}_x$ that are not shown are zero. We call $\mathbf{D}_x$ a differentiation matrix. Each row of $\mathbf{D}_x$ gives the weights of the 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. We are free to use whatever finite-difference formulas we like in each row. However, it makes sense to choose rows that are as similar as possible. Using second-order centered differences where possible and second-order one-sided formulas (see Table 5.4.2) at the boundary points leads to

The differentiation matrices so far are banded matrices, i.e., all the nonzero values are along diagonals close to the main diagonal.^[1]

10.3.2Second derivative¶

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

We have multiple choices again for $\mathbf{D}_{xx}$, and it need not be the square of any particular $\mathbf{D}_x$. As pointed out in Finite differences, squaring the first derivative is a valid approach but would place entries in $\mathbf{D}_{xx}$ farther from the diagonal than is necessary.

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.

10.3.3Spectral differentiation¶

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

1. Choose a node index set $S$ near node $i$.
2. Interpolate with a polynomial using the nodes in $S$.
3. 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:

then the resulting Chebyshev differentiation matrix has entries

where $c_0=c_n=2$ and $c_i=1$ for $i=1,\ldots,n-1$. Note that this matrix is dense. The simplest way to compute a second derivative is by squaring $\mathbf{D}_x$, as there is no longer any concern about the bandwidth 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).

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.4Exercises¶

1. (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''$?

2. (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.

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.

   (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.

4. ✍ Explain in detail how lines 23-24 in Function 10.3.2 correctly change the interval from $[-1,1]$ to $[a,b]$.

5. (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.

6. 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?