Figure A.2:

Gibbs energy contours for the pentane reactions as a function of the two reaction extents.

Code for Figure A.2

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
102
103
104
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

global K1 K2 P ext0 Gval dels unk zlast

deltag1 = -3.72; % kcal/mol
deltag2 = -4.49; % kcal/mol
T       = 400;   % K
P       = 2.5;   % atm
R       = 1.987; % cal/mol/K
K1 = exp(-deltag1*1000/(R*T));
K2 = exp(-deltag2*1000/(R*T));




% first optimize free energy
% min_z G(z) s.t. x_1,x_2 >0,  x_1 + x_2 < 0.5

Gval = 0;
w0 = [0.2; 0.2];
A = [1, 1];
b = 0.5;
[w, obj, info] = fmincon(@G, w0, A, b, [], [], [0;0], [0.5;0.5]);
w

if (info ~= 1)
  info, w, obj
end

% next prepare contour plot

Glevels = [0, -1, -2, -2.5, -2.5, -2.53, -2.55 -2.559];
nlev = length(Glevels);
alpha = [0.05, 0.05,   0.075, 0.05, 0.05, 0.05, 0.04, 0.02];
x0vec = [1e-3, 1e-3, 1e-3, .133, .133, .133, .133, .133; ...
          .02,  .15,  .35,   .3,  .3,   .3,    .3,   .3];
delx  = [0.01, 0.01, 0.01, 0.01, -0.01, 0, 0, 0; ...
            0     0,    0,    0,     0, 0, 0, 0];
sfinal = [9, 9, 9, 9, 9, 9, 9, 5];

for i = 1:nlev
  Gval = Glevels(i);
  %  find a point on the contour corresponding to s=0
  ext0 = x0vec(:,i);
  unk = [2];
  x0 = ext0(unk);
  [x, fval, info] = fsolve(@Gcontour, x0);
  fsolve_failed = info <= 0;
  infovec(i) = info;
  fcnvec(:,i) = Gcontour(x);
  ext = ext0;
  ext(unk) = x;
  zlast = ext;
  extinit(:,i) = ext;
  % continuation in arclength
  dels = alpha(i);
  svec = dels:dels:sfinal(i);
  narc =length(svec);
  z0 = ext + delx(:,i);
  unk = [1 2];
  for j = 1:narc
    s = svec(j);
    [z, fval, info] = fsolve('inflate',z0);
    fsolve_ok = info > 0;
    % stop when contour completed
    if (z(1) > 0 && z(2) > 0 && 0.5-(z(1)+z(2)) > 0 && fsolve_ok)
      ztable{i}(:,j) = z;
      % project initial guess along the arc for next point on contour
      z0 = 2*z - zlast;
      zlast = z;
    else
      break
    end
  end
  ztable(i) = {[ext, ztable{i}]'};
end
feasbound = [0, 0.5; 0.5, 0];  % store boundary of feasible region
save contour_pentane.dat ztable feasbound;


%plot results
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
  hold on;
  plot([0.5,0], [0, 0.5]); % one leg of feasible region
  for i = 1:nlev
    plot(ztable{i}(:,1), ztable{i}(:,2));
  end
% TITLE
end % PLOTTING

Gcontour.m


1
2
3
4
5
function retval = Gcontour(var)
  global K1 K2 P Gval ext0 unk
  w = ext0;
  w(unk) = var;
  retval = G(w);

G.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function retval = G(w)
  global K1 K2 P Gval
  x = w(1);
  y = w(2);
  z = 0.5 - x - y;
  p = 0.5 + z;
  if (x == 0)
    xlogx = 0;
  else
    xlogx = abs(x)*log(abs(x/p));
  end
  if (y == 0)
    ylogy = 0;
  else
    ylogy = abs(y)*log(abs(y/p));
  end
  if (z == 0)
    zlogz = 0;
  else
    zlogz = abs(z)*log(abs(z/p));
  end
  retval = -Gval -(x*log(K1) + y*log(K2)) + p*log(P) + ...
     2*zlogz + xlogx + ylogy;

inflate.m


1
2
3
function retval = inflate(z)
  global zlast dels
  retval = [ Gcontour(z) ; norm(z-zlast) - dels^2 ];