Stability region for TR-BDF2

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 the TR-BDF2 method.

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
In [3]:
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 [4]:
R = lambda z: (1 + 5/12.*z) / (1 - 7/12.*z + 1/12.*z**2)
plotS(R, axisbox=[-5,15,-10,10])

Boundary locus approach

Recall that for a LMM we can use the Boundary Locus method to find all points that could lie on the boundary of the stability region. We can try something similar for TR-BDF2. If a point $z$ lies on the boundary of $S$ then $|R(z)| = 1$ and so $R(z) = e^{i\theta}$ for some $\theta$ in $[0, 2\pi]$. For a LMM this gave an explicit formula for $z(\theta)$ that could be plotted in the complex plane. For TR-BDF2, since $R(z)$ is a rational function with a quadratic denominator, we can multiply by denominator and obtain a quadratic equation that would have to be solved for $z$ for each $\theta$. So there will be two curves $z_1(\theta)$ and $z_2(\theta)$ as $\theta$ varies.

The quadratic equation can be written as

$$ z^2 - (7 + 5\omega) + 12(1-\omega) = 0 $$

where $\omega = e^{-i\theta}$. Using the quadratic formula then gives the following code:

In [5]:
theta = linspace(-pi,pi,500)

omega = exp(-1j*theta)  # using 1j for the imaginary unit sqrt{-1}

# roots of quadratic:
s = sqrt((7+5*omega)**2 - 48*(1-omega))
z1 = ((7 + 5*omega) - s)/2
z2 = ((7 + 5*omega) + s)/2

plot(real(z1), imag(z1), 'r', label='root z1')
plot(real(z2), imag(z2), 'b--', label='root z2')

axisbox=[-5,15,-10,10]
xa, xb, ya, yb = axisbox
title('Boundary Locus for region of absolute stability')
grid(True)
plot([xa,xb],[0,0],'k')  # x-axis
plot([0,0],[ya,yb],'k')  # y-axis
legend()
axis('scaled')  # scale x and y same so that circles are circular
axis(axisbox);   # set limits

Note that the curves do cover the actual boundary of the stability region.

The vector of points theta above was chosen to go from $-\pi$ to $\pi$ so that each plot of $z_1$ and $z_2$ would result in a continuous curve. Note that when $\theta = 0$ the two roots are at $z_1 = 0$ and $z_2 = 12$, the midpoints of these curves. Hence if we instead choose theta to go from 0 to $2\pi$, the red curve, for example, would start at the origin, go move into the upper half plane, and then jump down to the complex conjugate of the uppermost red point, before returning smoothly to the origin. These jumps show up if we simply plot the curves and could be confusing, as seen below.

In [6]:
theta = linspace(0,2*pi,500)

omega = exp(-1j*theta)  # using 1j for the imaginary unit sqrt{-1}

# roots of quadratic:
s = sqrt((7+5*omega)**2 - 48*(1-omega))
z1 = ((7 + 5*omega) - s)/2
z2 = ((7 + 5*omega) + s)/2

plot(real(z1), imag(z1), 'r', label='root z1')
plot(real(z2), imag(z2), 'b--', label='root z2')

axisbox=[-5,15,-10,10]
xa, xb, ya, yb = axisbox
title('Boundary Locus for region of absolute stability')
grid(True)
plot([xa,xb],[0,0],'k')  # x-axis
plot([0,0],[ya,yb],'k')  # y-axis
legend()
axis('scaled')  # scale x and y same so that circles are circular
axis(axisbox);   # set limits
In [ ]: