Drawing stability domains in MATLAB

Grady Wright, 31st March 2015

Contents

Set up some variables for plotting.

LW = 'LineWidth'; lw = 1;
clr = [221 221 221]/255;

Stability domain for AB2

The main idea of tracing out stability domains is to apply the IVP solver, Adams-Bashforth (AB2)in this case, to $y'=\lambda y$. With a time-step of $k$ we get

$$y_{j+1} = y_{j} + \frac{3}{2} k\lambda y_{j} - \frac{1}{2} k\lambda y_{j-1} = y_{j} + \frac{\xi}{2}(3 y_j - y_{j-1}),$$

where we have set $\xi = k\lambda$. This is a linear recurrence equation for which the characteristic equation is

$$r^2 = r + \frac{\xi}{2}(3 r - 1).$$

The solution to the recurrence equation gives non-growing solutions when the roots of the characteristic equation are less than 1. To determine the region in the complex $\xi$-plane where this occurs, we first solve the above equation for $\xi$:

$$\xi = \frac{2r(r-1)}{3r-1}.$$

Next, we we plug in $r = \exp(i\theta)$ and let $\theta$ vary from $0$ to $2\pi$. This will trace out a curve in the complex $\xi$-plane indicating the boundary between growing (unstable) and non-growing (stable) solutions.

The code below illustrates the procedure and produces a plot of the stability domain. The grey-shaded region gives the values of $\xi$ where the scheme will be stable. Determining whether the stability region is inside or outside the closed curves requires picking a value of $\xi$ inside/outside and determining whether the resulting roots are less than one in modulus.

Define the unit circle in the complex plane

N = 1000;
th = linspace(0,2*pi,N);
r = exp(1i*th);

Solution of characteristic equation in terms of $\xi$

f = @(r) 2*r.*(r-1)./(3*r-1);

Evaluate $f$ at points on the unit circle and then plot the results:

xi = f(r);
plot(xi,'k-',LW,lw), hold on
fill(real(xi),imag(xi),clr)
plot([min(real(xi)) max(real(xi))],[0 0],'b--',LW,lw)
plot([0 0],[min(imag(xi)) max(imag(xi))],'b--',LW,lw)
axis tight, axis equal, hold off

Stability domains for other IVP methods are determined similarly. Below are several examples.

Stability domain for AB1 (Euler's method)

f = @(r) (r-1);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
fill(real(xi),imag(xi),clr)
plot([min(real(xi)) max(real(xi))],[0 0],'b--',LW,lw)
plot([0 0],[min(imag(xi)) max(imag(xi))],'b--',LW,lw)
axis tight, axis equal, hold off

Stability domain for AB3

f = @(r) 12*(r.^3-r.^2)./(23*r.^2 - 16*r + 5);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
fill(real(xi),imag(xi),clr)
plot([min(real(xi)) max(real(xi))],[0 0],'b--',LW,lw)
plot([0 0],[min(imag(xi)) max(imag(xi))],'b--',LW,lw)
axis tight, axis equal, hold off

Stability domain for AB4

Only the part in the left halfplane is included.

f = @(r) 24*(r.^4-r.^3)./(55*r.^3 - 59*r.^2 + 37*r - 9);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
fill(real(xi(real(xi)<=0)),imag(xi(real(xi)<=0)),clr)
plot([min(real(xi)) max(real(xi))],[0 0],'b--',LW,lw)
plot([0 0],[min(imag(xi)) max(imag(xi))],'b--',LW,lw)
axis tight, axis equal, hold off

Stability domain for BDF2

For all values outside of the circle the scheme is stable.

f = @(r) (3*r.^2 - 4*r + 1)./(2*r.^2);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
set(gca,'Color',clr);
fill(real(xi),imag(xi),'w')
axis([-4 4 -3 3])
axis equal
ax = axis;
plot([ax(1) ax(2)],[0 0],'b--',LW,lw)
plot([0 0],[ax(3) ax(4)],'b--',LW,lw)
hold off

Stability domain for BDF3

For all values outside of the circle the scheme is stable.

f = @(r) (11*r.^3 - 18*r.^2 + 9*r - 2)./(6*r.^3);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
set(gca,'Color',clr);
fill(real(xi),imag(xi),'w')
axis([-4 7 -4 4])
axis equal
ax = axis;
plot([ax(1) ax(2)],[0 0],'b--',LW,lw)
plot([0 0],[ax(3) ax(4)],'b--',LW,lw)
hold off

Stability domain for BDF4

For all values outside of the circle the scheme is stable.

f = @(r) (25*r.^4 - 48*r.^3 + 36*r.^2 - 16*r + 3)./(12*r.^4);
xi = f(r);
plot(xi,'k-',LW,lw), hold on
set(gca,'Color',clr);
fill(real(xi),imag(xi),'w')
axis([-4 12 -7 7])
axis equal
ax = axis;
plot([ax(1) ax(2)],[0 0],'b--',LW,lw)
plot([0 0],[ax(3) ax(4)],'b--',LW,lw)
hold off