## Code for Figure 2.26

### 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% 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);

nu2 = (z.^2 - z)./((3*z - 1)/2);

nu3 = (z.^3 - z.^2)./((5-16*z + 23*z.^2)/12);

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
plot(real(nu2),imag(nu2),'r-','LineWidth',2, ...
real(nu3),imag(nu3),'b-','LineWidth',2, ...
real(nu4),imag(nu4),'g-','LineWidth',2, ...
[xa xb],[0 0],'k-','LineWidth',2, ...
[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