2. Local interpolation using nodal values : finite differences

Ronojoy Adhikari (rjoy@imsc.res.in) http://rjoy.in/group | The latest version of this article is available at https://github.com/ronojoy/pyddx-lectures

Polynomial interpolants approximate functions by their values at grid points. If the function values are known at the grid points $x_0, x_2, \ldots x_N$ the Lagrange interpolating polynomial is

$$ p_{N}(x) = \sum_{k=0}^N f(x_k) F_k(x) $$

where

$$ F_k(x) = \frac{\prod_{j=0 \atop j\neq k}^N(x-x_j)}{\prod_{j=0 \atop j\neq k}^N(x_k - x_j)}, \\,\\, k = 0, 1, \ldots N $$

are $N$-th degree polynomials satisfying $F_k(x_j) = \delta_{kj}$. By construction, the interpolating polynomial attains the function values at the nodes, $p_N(x_k) = f(x_k)$, but its behaviour between the nodes depends on that of $F_k(x)$. Derivatives of the interpolant, in particular, are determined by derivatives of the $F_k(x)$,

\begin{align} \frac{d^m}{dx^m}p_N(x) = \sum_{k=0}^N \frac{d^m}{dx^m} f(x_k) F_k(x). \end{align}

The polynomial approximations for the derivatives of $f(x)$ at $x = \xi$ using the nodal values $f(x_k)$ are obtained by evaluating the derivatives of the interpolant at $ x = \xi$,

\begin{align} f^{(m)}(\xi) = \frac{d^m}{dx^m}p_N(x)|_{\xi} = \sum_j^N \left[ \frac{d^m}{dx^m} F_k(x) \right]_{\xi} \cdot f(x_k). \end{align}

Explain here Fornberg's algorithm for calculating the weights.

Then, specialise to the case where the points of evaluation are identical to the points where the function values are known, and then say that it can be expressed as differentiation matrices.

\begin{align} f^{(m)}(x_j) = \frac{d^m}{dx^m}p_N(x)|_{x_j} = \sum_j^N \left[ \frac{d^m}{dx^m} F_k(x) \right]_{x_j} \cdot f(x_k) \end{align}

and the operation can be thought of as a differentiation matrix ${\bf D}^{(m)}$ multiplying the vector of function values ${\bf f}$ to give a vector of derivative values ${\bf f}^{(m)}$,

$$ f^{(m)}_j = f^{(m)}(x_j) = \sum_{k=0}^N D^{(m)}_{jk}f(x_k) $$

where elements of the differentiation matrix are

$$ D^{(m)}_{jk} = \left[ \frac{d^m}{dx^m} F_k(x) \right]_{x_j}. $$

Introduction

The finite-difference method for solving partial differential equations rests on low-order polynomial interpolation of the function, typically over a regular spaced set of "collocation" points. Expressions for the function gradient are obtained by differentiating the interpolant and evaluating the derivative at the collocation points. To each choice of interpolant corresponds a finite-difference formula for the gradients. Finite-difference gradients inherit the error between the interpolant and the function and this error is equal to the order of the polynomial interpolant.

To make these notions definite, consider the function $y(x)$ whose values are given at equally spaced "collocation" points $x_1, x_2, \ldots x_N$ with $x_{i+1} - x_i = h$ for every $i$. An interpolant, $p(x)$, is a function that coincides with $y(x)$ at the collocation points, that is

$$ p(x_i) = y(x_i) \equiv y_i \quad \forall i $$

but may differ from $y(x)$ at other points in the interval $(x_1, x_N)$. Local polynomial interpolation constructs the interpolant from an union of local polynomial. A $m$-th order polynomial contains $m+1$ unknown coefficients which require that many conditions to be determined. The conditions are usually taken to be that the interpolant coincides with the function at the collocation points. For example, the trhee unknown coefficients of a local quadratic interpolant $p_i(x)$ are determined by requiring

$$ p_i(x_i) = y_i, \quad p_i(x_{i-1}) = y_{i-1}, \quad p_i(x_{i+1}) = y_{i+1} $$

The solution for $p_i(x)$ is provided by the famous Lagrange interpolation formula

$$ p_i(x) = y_{i-1}\frac{x - x_{i}}{x_{i-1} - x_{i}}\frac{x - x_{i+1}}{x_{i-1} - x_{i+1}} + y_{i}\frac{x - x_{i-1}}{x_{i} - x_{i-1}}\frac{x - x_{i+1}}{x_{i} - x_{i+1}} + y_{i+1}\frac{x - x_{i-1}}{x_{i+1} - x_{i-1}}\frac{x - x_{i}}{x_{i+1} - x_{i}} $$

The finite-difference derivative of $y(x)$ is now simply the derivative of the interpolant, evaluated at the collocation point :

$$ y^{\prime}(x_i) = p_i^{\prime}(x) |_{x_i} $$

This yields the well-known central difference formula for the first derivative

$$ y^{\prime}(x_i) = \frac{y_{i+1} - y_{i-1}}{ h^2 }$$

Extending this recipe to second derivatives gives the well-known central difference formula for the Laplacian,

$$ y^{\prime\prime}(x_i) = \frac{y_{i+1} -2 y_i + y_{i-1}}{ h^2 }$$

Both the above formulas are $O(h^2)$ accurate as they are obtained by differentiating a quadratic interpolant.

See the notes by Detlev Maurer for more general interpolation : https://people.fh-landshut.de/~maurer/numeth/numeth.html

Multidimensional Hermite polynomials : http://www.jstor.org/discover/10.2307/2004829?uid=3739808&uid=2&uid=4&uid=3739256&sid=21102813839791

Sampling theorem

A central result in interpolation is the Nyquist-Shannon sampling theorem, which Shannon states as

If a function x(t) contains no frequencies higher than B hertz, it is completely determined by giving its ordinates at a series of points spaced 1/(2B) seconds apart.

In other words, if a function is sufficiently smooth, it can be reconstructed completely from a finite number of its samples - a remarkable result!

Polynomial interpolation

Finite difference methods follow from the solution of the following interpolation problem. Given a set of grid points $x_0, x_1, \ldots, x_n$ and a point $ x = \xi$, not necessarily coincident with the grid points, find the weights $c_{i, j}^k$ such that the approximations

$$ \frac{d^pf}{dx^p}|_{x=\xi} \approx \sum