Figure 2.26:

Stability regions for Adams-Bashforth methods; \dot {x}=\lambda x.

Code for Figure 2.26

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
% Program to plot stability plots for Adam Bashforth methods
% For details on plotting stability region for Linear Multistep Method
% (LMM), see pg 162 of Randall J. LeVeque book on " Finite Difference
% Methods for Ordinary and Partial Differential Equations

% The characteristic polynomial of a LMM is pi(s,nu)=rho(s)-nu*sigma(s)
% where nu = lambda*dt
% A linear multistep formula is absolutely stable for a particular value of
% nu if and only if all the zeros of the stability polynomial pi satisfy
% |s| < 1 and any zero with |s| = 1 is a simple root.

% We have to find the value of nu s.t. s = exp(i*theta) is a root of
% stability polynomial pi.

% We demonstrate here stability curve for AB2, AB3 and AB4 methods.

% Setting grid points

axisbox = [-1.5 1.5 -1 1];
xa = axisbox(1); xb = axisbox(2); ya = axisbox(3); yb = axisbox(4);
npts  = 50;
theta = linspace(0, 2*pi, 2*npts+1);
z=exp(i*theta);

% 2nd order Adam Bashforth
nu2 = (z.^2 - z)./((3*z - 1)/2);

% 3rd order Adam Bashforth
nu3 = (z.^3 - z.^2)./((5-16*z + 23*z.^2)/12);

% 4th order Adam Bashforth
nu4 = (z.^4 -z.^3)./((55*z.^3 -59*z.^2 +37*z -9)/24);

% chop off intersecting loops for 4th order AB
  for k = 1: length(nu4)-1
    z = [real(nu4); imag(nu4)];
    iloop = 0;
    for j = k+2 : length(nu4)-1
      % find first place where z(k) intersects the rest of the path
      lam = inv([z(:,k)-z(:,k+1), z(:,j+1)-z(:,j)])*(z(:,j+1)-z(:,k+1));
      if (lam >=0 & lam <=1)
	iloop = k+1:j;
	zint = lam(1)*z(:,k)+(1-lam(1))*z(:,k+1);
	break
      end
    end
    if (iloop ~= 0)
      % chop out nu4(iloop) and replace with single interpolated value zint
      zcp = complex(zint(1),zint(2));
      nu4(iloop(1)) = zcp;
      iloop(1) = [];
      nu4(iloop) = [];
    end
  end

figure()
clf
hold on
plot(real(nu2),imag(nu2),'r-','LineWidth',2)
plot(real(nu3),imag(nu3),'b-','LineWidth',2)
plot(real(nu4),imag(nu4),'g-','LineWidth',2)
plot([xa xb],[0 0],'k-','LineWidth',2)
plot([0 0],[ya yb],'k-','LineWidth',2)
legend('AB2','AB3','AB4')
set(gca,'FontSize',15)
title('Region of absolute stability of AB methods','FontSize',15)

nus = [real(nu2); imag(nu2); real(nu3); imag(nu3)]';
nus4 = [real(nu4);imag(nu4)]';
figure()
plot(z(1,:),z(2,:), '-o')

figure()
plot(real(nu4), imag(nu4), '-o')

save "AB_stability.dat"  nus nus4