Stability regions of some 1-step methods

AMath 586, Spring Quarter 2019 at the University of Washington. For other notebooks, see Index.html or the Index of all notebooks on Github.

Plot the region of absolute stability for various 1-step methods.

The general approach is to apply the method to $u' = \lambda u$ with time step $k$ to obtain $U^{n+1} = R(z) U^n$, where $R(z)$ is a rational function of $z=k\lambda$ (a polynomial for an explicit method). Then evaluate $|R(z)|$ on a grid of points in the complex plane and do a filled contour plot that shows the regions where $|R(z)| \leq 1$.

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
seterr(divide='ignore', invalid='ignore')  # suppress divide by zero warnings
Out[2]:
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
In [3]:
from ipywidgets import interact
In [4]:
def plotS(R, axisbox = [-10, 10, -10, 10], npts=500):  
    """
    Compute |R(z)| over a fine grid on the region specified by axisbox
    and do a contour plot with contourf (filled contours)
    to show the region of absolute stability.
    """

    xa, xb, ya, yb = axisbox
    x = linspace(xa,xb,npts)
    y = linspace(ya,yb,npts)
    X,Y = meshgrid(x,y)
    Z = X + 1j*Y
    Rval = R(Z)
    Rabs = abs(Rval)
    
    # plot interior, exterior, as green and white:
    levels = [-1e9,1,1e9]
    CS1 = contourf(X, Y, Rabs, levels, colors = ('g', 'w'))

    # plot boundary as a black curve:
    CS2 = contour(X, Y, Rabs, [1,], colors = ('k',), linewidths = (2,))
    
    title('Region of absolute stability')
    grid(True)
    plot([xa,xb],[0,0],'k')  # x-axis
    plot([0,0],[ya,yb],'k')  # y-axis

    axis('scaled')  # scale x and y same so that circles are circular
    axis(axisbox)   # set limits
In [5]:
R = lambda z: 1+z
plotS(R, axisbox=[-5,5,-5,5])

Theta Method

The Theta Method is

$$U^{n+1} = U^n + k[(1-\theta)f(U^n) + \theta f(U^{n+1}]. $$

The right hand side is a convex combination of $f(U^n)$ and $f(U^{n+1})$. Special cases are:

  • Forward Euler: $\theta = 0$,
  • Trapezoid: $\theta = 1/2$,
  • Backward Euler: $\theta = 1$.

Note how the region of absolute stability varies with $\theta$, and that the method is A-stable for $\theta \geq 1/2$.

In [6]:
def plotS_theta(theta):
    R = lambda z: (1. + (1-theta)*z) / (1-theta*z)
    plotS(R, npts=50)  # use fewer points so interact works well
In [7]:
interact(plotS_theta, theta=(0,1,.1));

Taylor series methods

When applied to $u' = \lambda u$, the $p$th order Taylor series method gives $R(z)$ equal to the first $p+1$ terms of the Taylor series expansion of $e^z = \sum_{j=0}^\infty z^j/j!$.

In [8]:
def plotS_TS(r):
    def R(z):
        # return Rz = 1 + z + 0.5z^2 + ... + (1/r!) z^r
        Rz = 1.
        term = 1.
        for j in range(1,r+1):
            term = term * z/float(j)
            Rz = Rz + term
        return Rz
            
    plotS(R, npts=500)  # use fewer points so interact works well
    title('Taylor series method r = %i' % r)
In [9]:
interact(plotS_TS, r=(1,20,1));