Stability region for a simple Predictor-Corrector method

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 Predictor-Corrector method based on Forward and Backward Euler,

\begin{align*} &\hat U^0 = U^n + kf(U^n)\\ &\text{for} ~j = 0,~1,~\ldots,~N-1\\ &\qquad \hat U^{j+1} = U^n + kf(\hat U^j)\\ &\qquad \text{end}\\ &U^{n+1} = \hat U^N. \end{align*}

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$. This gives

$$ R(z) = 1 + z + z^2 + \cdots + z^{N+1} = \frac{1-z^{N+2}}{1-z}. $$

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]:
from ipywidgets import interact
In [4]:
def plotS(R, axisbox = [-2, 2, -2, 2], 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]:
def plotS_PC(N):
    def R(z):
        # return Rz = 1 + z + z^2 + ... + z^(N+1) = (1 - z^(N+2)) / (1-z)
        return (1 - z**(N+2)) / (1-z)
     
    plotS(R, npts=1000)

    title('PC method N = %i' % N)
In [6]:
figure(figsize=(14,4))
subplot(1,4,1)
plotS_PC(5)
subplot(1,4,2)
plotS_PC(10)
subplot(1,4,3)
plotS_PC(20)
subplot(1,4,4)
plotS_PC(50)

The limiting region

As $N \rightarrow \infty$, the stability region appears to approach the intersection of the unit circle and the stability region of the backward Euler method, which is the exterior of the circle of radius 1 centered at $z=1$. This makes sense since (a) the P-C iteration converges only if $|z| \leq 1$ and (b) if it does converge then it converges to the value $U^{n+1}$ that solves the implicit equation defining the backward Euler method.

The blobs on the unit circle

Note that the stability region seems to also include an increasing number of small blobs equally spaced around the unit circle. This is because, for any finite $N$, the function $R(z)$ is given by

$$ R(z) = \frac{1-z^{N+2}}{1-z}. $$

If $|z|<1$ then at points outside the stability region of backward Euler the limiting value $1/(1-z)$ is greater than 1 in modulus, but for finite $N$ the quantity $1 - z^{N+2}$ in the numerator vanishes if $z^{N+2} = 1$, i.e., at the $N+2$ roots of unity $\exp({2\pi ij/(N+2)})$ for $j = 1, 2, \ldots,~ N+2$. These points are equally spaced around the unit circle and the plot below confirms that the blobs are centered around these points, which are expressed as $\omega^j$ where $\omega = \exp({2\pi i/(N+2)})$.

In [7]:
figure(figsize=(10,10))

N = 50
plotS_PC(N)

omega = exp(2*pi*1j/(N+2))
roots_of_unity = omega**range(1,N+3)
plot(real(roots_of_unity), imag(roots_of_unity), 'rx', 
     markersize=5, label='roots of unity')
legend(fontsize=15)
Out[7]:
<matplotlib.legend.Legend at 0x1270fad90>

Note also that half way between each of the roots of unity, at $z = \exp({2\pi i(j+1/2)/(N+2)})$, we have $z^{N+2} = \exp(\pi i) = -1$ and so $R(z) = 2/(1-z)$. Near these points $|R(z)| >1$ which explains the scalloped shape of the stability region where it extends out beyond the unit circle a little bit around the roots of unity in the plot above, but then inside the unit circle in between.

These blobs vanish in the limit as $N \rightarrow \infty$.

Note that letting $N\rightarrow \infty$ does not entirely make sense in the context of the predictor-corrector method, since this means taking an infinite number of iterations of the corrector in each time step. In practice we would iterate until the fixed point iteration converges to within some tolerance. This can only happen if $|z|<1$ so the stability region for the full method is contained in the open unit disc.