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.

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

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

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

to be applied to a vector of function values at the nodes (10.3.1).

10.3.2Second derivative¶

In a TPBVP, we will need to take the second derivative of the unknown solution. One option is to apply a first derivative twice, that is, to multiply the value vector $\mathbf{f}$ on the left by $\mathbf{D}_x^2$. This is usually not the best option, however (see Finite differences). Instead, the following is usually a better choice:

10.3.3Implementation¶

Together, the matrices (10.3.6) and (10.3.10) 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.4Spectral 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 the nodes each time.

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

Function 10.3.2 returns the matrices from Definition 10.3.4. This 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.12) (see Exercise 10.3.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.5Exercises¶