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.
In [ ]:
import sys, os
# pyddx is installed in ~/Numerics/pyddx
sys.path.append(os.path.expanduser('~/Numerics'))
# import modules from pyddx
from pyddx.fd import fdsuite as fds
from pyddx.sc import dmsuite as dms
In [22]:
NN = 2**arange(3,14)
fig, (ax1, ax2) = subplots(1, 2, figsize=(12, 4))
for N in NN:
h = 2.*pi/N; x = -pi + arange(1,N+1).T*h
u = exp(sin(x)); uprime = cos(x)*u
D3pt = fds.cdif(N, 1, 3, h); D5pt = fds.cdif(N, 1, 5, h)
error3pt = linalg.norm(D3pt.dot(u) - uprime,inf); error5pt = linalg.norm(D5pt.dot(u) - uprime, inf);
#print "3pt error =",error3pt, "\t ", "5pt error =", error5pt, " ", "N =", N
ax1.loglog(N,error3pt,'.',markersize=15); ax2.loglog(N,error5pt, '*', markersize=15)
# plot both sets of results
ax1.loglog(NN, NN**(-2.), '--')
ax2.loglog(NN, NN**(-4.), '--')
ax1.text(105.,2.e-7,'N^-2',fontsize=18)
ax2.text(105.,2.e-14,'N^-4',fontsize=18)
ax1.grid(); ax2.grid();
ax1.set_xlabel('N'), ax1.set_ylabel('Error')
ax2.set_xlabel('N'), ax2.set_ylabel('Error')
ax1.set_title('3 point stencil accuracy')
ax2.set_title('5 point stencil accuracy')
Out[22]:
In [ ]: