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
These lectures will introduce you to pyddx, a differentation matrix suite in Python, and its use for solving ordinary and partial differential equations.
Differentiation matrices provide accurate and easy to implement methods for solving differential equations. Differentiation matrices are derivatives of interpolants, functions chosen to approximate the desired solution of a differential equation. Each interpolation method, that is, a choice of interpolation function and a choice of the manner in which the solution is represented by a finite number of values, yields an approximation for the solution.
The pyddx library contains differentiation matrices that follow from three kinds of interpolation schemes : (a) local low order polynomial interpolation (b) local low order interpolation using isotropic polynomials and (c) global interpolation using orthogonal functions. These yield, respectively, (a) standard finite difference methods, (b) kinetic finite difference methods and (c) spectral collocation methods. Taken together, these differentiation matrices provide flexible, easily implemented, and accurate methods to solve partial differentiation equations in simple domains.
Interpolation is a key idea that permeates much of numerical analysis and forms the basis of numerical approximation schemes for derivatives and integrals. The basic idea is to approximate a function $f(x)$, whose values $f_j = f(x_j)$, are known at a finite number of collocation points $x_1, x_2, \ldots, x_N$. This is accomplished through the interpolant $p(x)$
$$ f(x) \approx p(x) = \sum_{j=1}^N \frac{\alpha(x)}{\alpha(x_j)}\phi_j(x) f_j $$where $\alpha(x)$ are weight functions and $\phi_j(x)$ are basis functions of the interpolation scheme. For $p(x)$ to be an interpolant, it is necessary that it agrees with the values of the function at the collocation points, $ p(x_j) = f_j$, which requires the basis functions to satisfy $\phi_j(x_k) = \delta_{jk}$. Common choices of interpolants include piecewise polynomials of low order and orthogonal polynomials. We will see below how each choice of collocation points $x_j$, weights $\alpha(x)$ and basis functions $\phi_j(x)$ give rise to different approximations for derivatives.
Given the interpolant, derivatives of $f(x)$ can be computed by taking derivatives of $p(x)$ and evaluating them at the collocation points,
$$ f^{(n)}(x_k) \approx \sum_{j=1}^N \frac{d^n}{dx^n}\left[ \frac{\alpha(x)}{\alpha(x_j)}\phi_j(x)\right]_{x=x_k} f_j $$Since both interpolation and differentiation are linear operations they can be represented by matrices acting on vectors. Thus the differentiation above can be expressed as multiplication by a differentiation matrix $\mathbf{D}^{(l)}=D^{(l)}_{kj}$ acting on the function vector $\mathbf{f} = (f_1, f_2, \ldots, f_N)$ to produce the $n$-th derivative vector $\mathbf{f}^n = (f^{n}(x_1), f^{n}(x_2), \ldots, f^{n}(x_n))$,
$$ \mathbf{f}^n = \mathbf{D}^{(l)}\cdot \mathbf{f} $$where the differentiation matrix is given by
$$ \mathbf{D}^{(l)}=D^{(l)}_{kj} = \frac{d^n}{dx^n}\left[ \frac{\alpha(x)}{\alpha(x_j)}\phi_j(x)\right]_{x=x_k} $$This work benefitted from the contributions of many previous authors. In particular, the monographs by Nick Trefethen (Spectral Methods in Matlab) and Bengt Fornberg (A practical guide to pseudospectral methods) have been very helpful. Weidemann and Reddy's differentiation matrix suite and Nick Trefethen's spectral collocation codes, both in Matlab, have been sources of inspiration. Robert Piche's short and clear examples on solving differential equations using differentiation matrices have been incorporated in the accompanying lecture notes. Several authors have attempted to provide spectral collocation matrices in Python, for example, chebpy and another-chebpy - pyddx incorporates some of their code for spectral differentiation matrices. Needless to say, none of this would have been possible without the collective effort of the community that has given us Numpy, Scipy, Ipython and other wonderful software tools.