Order stars

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 relative stability, also called the order star, 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)/e^z|$ on a grid of points in the complex plane and do a filled contour plot that shows the regions where $|R(z)/e^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 plotOS(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 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)
    Rrel = abs(Rval / exp(Z))
    
    # plot interior, exterior, as green and white:
    levels = [-1e9,1,1e9]
    CS1 = contourf(X, Y, Rrel, levels, colors = ('g', 'w'))

    # plot boundary as a black curve:
    CS2 = contour(X, Y, Rrel, [1,], colors = ('k',), linewidths = (2,))
    
    title('Order Star')
    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
plotOS(R, axisbox=[-5,5,-5,5])

Theta method

In [6]:
def plotOS_theta(theta):
    R = lambda z: (1. + (1-theta)*z) / (1-theta*z)
    figure(figsize=(6,6))
    plotOS(R, npts=200)  # use fewer points so interact works well
    title("Order star for theta-method with theta = %4.2f" % theta)

For $\theta=1/2$ this is the Trapezoid method, which is second order accurate, and so the order star has 3 sectors inside the region of relative stability near $z=0$:

In [7]:
plotOS_theta(0.5)

For all other $\theta$, the method is only first order accurate and there are only 2 sectors inside/outside the order star near $z=0$:

In [8]:
interact(plotOS_theta, theta=(0,1,.1));

Taylor series methods

For this class of methods we can easily increase the order and observe how the structure near the origin varies:

In [9]:
def plotOS_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
            
    figure(figsize=(6,6))
    plotOS(R, npts=300)
    title('Taylor series method r = %i' % r)
In [10]:
interact(plotOS_TS, r=(1,20,1));

Note that the way this computation is done it is subject to a lot of cancellation near $z=0$, so plots above do not look correct in this region, particularly for $r>15$ or so -- there should be $r+1$ green sectors and $r+1$ white sectors approaching the origin.

Here's a better version for the particular case of high-degree Taylor series methods:

In [11]:
def plotOS_TS(r):
    def R(z):
        # return Rz = 1 + z + 0.5z^2 + ... + (1/r!) z^r
        from math import factorial
            
        Rz = 1.
        term = 1.
        for j in range(1,r+1):
            term = term * z/float(j)
            Rz = Rz + term
        
        # for small z when r is large, 
        # it's better to compute the remainder
        remainder = 0.
        term = 1./factorial(r) 
        for j in range(r+1, 2*r):
            term = term * z/float(j)
            remainder = remainder + term
        remainder *= z**r 
        
        # Define this so that contours of |R(z)/exp(z)| = 1 
        # look right near origin:
        Rz_smallz = where(real(remainder/exp(z))>0., 0.5*exp(z), 2*exp(z))           
        Rz = where(abs(z**(r+1)/factorial(r+1)) > 1e-10, Rz, Rz_smallz)
        
        return Rz
            
    figure(figsize=(6,6))
    plotOS(R, npts=500)
    title('Taylor series method r = %i' % r)
In [12]:
interact(plotOS_TS, r=(1,30,1));