Figure 2.27:

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

Code for Figure 2.27

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
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
% Program to plot stability plots for predictor-corrector methods

% A linear multistep formula is: x+ = G(w) x  with w = a*Delta t
% for differential equation xdot = ax with stepsize Delta t
%
% The eigenvalues of G determine the stability of the method;
% The method is absolutely stable iff max (abs( eig G(w))) < 1

% We find the value of w s.t. z = exp(i*theta) is an eigenvalue
% of G for theta in [0, 2*pi].  Actually, [0, 4*pi] to cover entire
% curve in w.

% We demonstrate here stability curve for APC2, APC3 and APC4 methods.
% NOTE: APCn uses an n-1st order predictor and nth order corrector
%

% Setting grid points

% pass parameters to function, jbr, 6/28/2018

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

% predictor coefficients, 1st, 2nd, 3rd, 4th order
p{1} = [1, 0, 0, 0];
p{2} = [3/2, -1/2, 0, 0];
p{3} = [23/12, -4/3, 5/12, 0];
p{4} = 1/24*[55, -59, 37, -9];
% corrector coefficients, 1st, 2nd, 3rd, 4th order
c = cell(4);
c{1} = [1, 0, 0, 0];
c{2} = [1/2, 1/2, 0, 0];
c{3} = 1/12*[5, 8, -1, 0];
c{4} = 1/24*[9, 19, -5, 1];

par = struct();
par.p = p; par.c = c;


PCsvec = [2,3,4];
npcs = length(PCsvec);
wsol(1:npcs) = {zeros(npts,1)};
%optimset("ComplexEqn", "on");
for j = 1:npcs
  PCs = PCsvec(j);
  par.PCs = PCs;
  w0 = 0 + 0*i;
  for k = 1:length(theta)
    lam = lamvec(k);
    par.lam = lam;
    [w,fval,info] = fsolve(@(w) characteristic(w, par), w0);
    if (info <= 0 )
      j, k, info, w, w0
    end
    w0 = w;
    wsol{j}(k) = w;
  end
end

% check for intersecting loops
for ind = 1:npcs
  for k = 1: length(wsol{ind})-1
    z = [real(wsol{ind})'; imag(wsol{ind})'];
    iloop = 0;
    for j = k+2 : length(wsol{ind})-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 wsol(iloop) and replace with single interpolated value zint
      zcp = complex(zint(1),zint(2));
      wsol{ind}(iloop(1)) = zcp;
      iloop(1) = [];
      wsol{ind}(iloop) = [];
    end
  end
end

figure()
hold on
for j = 1:npcs
  plot(real(wsol{j}), imag(wsol{j}),'-o')
end
%axis (axisbox)

% rewrite real and imaginary parts for plotting
store = cell(npcs,1);
for k = 1: npcs
  tmp = cell2mat(wsol(k));
  store(k) = {[real(tmp), imag(tmp)]};
end
save ("APC_stability.dat", "store")

characteristic.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
function retval = characteristic(w, par)
  % note, the predictor order is one less than corrector
  cc = par.c{par.PCs}; pp=par.p{par.PCs-1};
      a(1) = 1 + w*( cc(1)*(1+w*pp(1)) + cc(2) );
      a(2) = w*( cc(1)*w*pp(2) + cc(3) );
      a(3) = w*( cc(1)*w*pp(3) + cc(4) );
      a(4) = w*cc(1)*w*pp(4);
  order = length(a);
  G = [a; eye(order-1), zeros(order-1,1)]-par.lam*eye(order);
  retval = det(G);