CBE 504 HW1

hw1.m

N = 10;
k = 1;
P0 = [zeros([N, 1]) ; 1];
A = irreversible(N);
f = @(t,P)mastereqn(P,A);
[tvec,P]=ode45(f, [0,25], P0);
tvecsize = size(tvec);
num_timesteps = tvecsize(1);

%  % Uncomment to check probability sums:
%  for i=[1:num_timesteps]
%      sum = 0;
%      for j=[1:N+1]
%          sum += P(i,j);
%      end
%      sum % this should print a bunch of varieties of 1.
%  end

t = repmat(tvec, 1, N+1);
tmin_plot = 0;
tmax_plot = 4;
nvec = [0:N];
n = repmat(nvec, num_timesteps, 1);

fig1 = figure(1);
%replace surf() with stem3() to get a more obviously discrete 3D plot (more like a bar chart)
stem3(n, t, P);
axis([0 N tmin_plot tmax_plot 0 1]);
xlabel('n (number of surviving molecules)');
ylabel('t [seconds]');
zlabel('P(n,t)');
title(strcat('Time-dependent probability distribution for N=',num2str(N),' molecules undergoing irreversible decomposition.'));
outsideaxes = axes('Position',[-.9 -.9 1 1],'Visible','off','YColor','w');
hold on;
set(gcf,'CurrentAxes', outsideaxes);
text(0.0001,0.0001,'CBE 504, HW1-1a, Tom Bertalan')
print(fig1,'hw1-1a.png');

% Use the resulting numerical solution to calculate the time-dependent mean and standard deviation of N in the system. Compare these numerical results with analytical results derived in class.

means = zeros(num_timesteps,1);
sdevs = zeros(num_timesteps,1);
for i = [1:num_timesteps]
    sum_of_X = 0;
    for j = [1:N+1]
        means(i) = means(i) + P(i,j) * nvec(j);
    end
    sum_of_sqdiffs = 0;
    for j = [1:N+1]
%          fprintf('sum_of_sqdiffs = %f + ( %f * %f - %f )^2;\n', sum_of_sqdiffs, nvec(j), P(i,j), means(i) );
        sum_of_sqdiffs = sum_of_sqdiffs + P(i,j) * ( nvec(j) - means(i) )^2;
    end
    sdevs(i,1) = sum_of_sqdiffs;
end

expectations = zeros(num_timesteps,1);
variances = zeros(num_timesteps,1);
for i = [1:num_timesteps]
    expectations(i) = N*exp(-1*k*tvec(i,1));
    variances(i)    = expectations(i) * (1 - exp(-1*k*tvec(i,1)) );
end


fig2 = figure(2);
lax = axes('YAxisLocation','left',  'YColor','k');
hold on; % hold on needs to be done after each axis creation ... 
rax = axes('YAxisLocation','right', 'YColor','b');
hold on; % ... not just per figure.
explot   =    plot(lax, tvec,expectations,'0');
mscatter = scatter(lax, tvec,means,[],[0 0 0]);
varplot  =    plot(rax, tvec,variances,'b');
sscatter = scatter(rax, tvec,sdevs,[],[0 0 1]);
legend(lax, 'expectation (curve is analytical; dots are data mean)', 'location', 'NorthEast');
legend(rax, 'variance (curve is analytical; dots are data variance)', 'location', 'East');
xlabel('t [sec]')
ylabel(lax, 'expectation');
ylabel(rax, 'variance');
title(strcat('Analytical and Numerical variance and expectation for number of survivors among N=',num2str(N),' molecules undergoing irreversible decomposition.'));
outsideaxes = axes('Position',[-.9 -.9 1 1],'Visible','off','YColor','w');
hold on;
set(gcf,'CurrentAxes', outsideaxes);
text(0.0001,0.0001,'CBE 504, HW1-1b, Tom Bertalan')
print(fig2,'hw1-1b.png')

%  Number 2 - using Matrix Exponentiation instead of ODE45 (Runge-Kutta, probably).
N=30;
k1=1;
k2=3;
A = reversible(N,[k1,k2])
t0=0;
dt=.1;
tf=10;
tvec = transpose([t0:dt:tf]);
tvecsize=size(tvec);
num_timesteps = tvecsize(1);
P = zeros([num_timesteps,N+1]);
P(1,N+1)=1;
for i=[2:num_timesteps]
    P(i,:) = transpose(mx_exp_solve(A,tvec(i),transpose(P(1,:))));
end

%  % Uncomment to check probability sums:
%  for i=[1:num_timesteps]
%      sum = 0;
%      for j=[1:N+1]
%          sum += P(i,j);
%      end
%      sum % this should print a bunch of varieties of 1.
%  end

t = repmat(tvec, 1, N+1);
tmin_plot = 0;
tmax_plot = 4;
nvec = [0:N];
n = repmat(nvec, num_timesteps, 1);

fig3 = figure(3);
%replace surf() with stem3() to get a more obviously discrete 3D plot (more like a bar chart)
stem3(n, t, P);
axis([0 N tmin_plot tmax_plot 0 1]);
xlabel('n (number of molecules in state 1)');
ylabel('t [seconds]');
zlabel('P(n,t)');
title(strcat('Time-dependent probability distribution for N=',num2str(N),' molecules undergoing reversible isomerization.'));
outsideaxes = axes('Position',[-.9 -.9 1 1],'Visible','off','YColor','w');
hold on;
set(gcf,'CurrentAxes', outsideaxes);
text(0.0001,0.0001,'CBE 504, HW1-2a, Tom Bertalan')
print(fig3,'hw1-2a.png');

tvec = [0; 0.1/(k1+k2); 1/(k1+k2); 2/(k1+k2); 10/(k1+k2) ; 100000/(k1+k2) ];
tvecsize=size(tvec);
num_timesteps = tvecsize(1);
P = zeros([num_timesteps,N+1]);
P(1,N+1)=1;
for i=[2:num_timesteps]
    P(i,:) = transpose(mx_exp_solve(A,tvec(i),transpose(P(1,:))));
end
nvec = [0:N];


fig4 = figure(4);
ax = axes('YAxisLocation','left',  'YColor','k');
hold on;
plot(ax, nvec, P(1,:), '-+r')
plot(ax, nvec, P(2,:), '-ob')
plot(ax, nvec, P(3,:), '-dg')
plot(ax, nvec, P(4,:), '-sy')
plot(ax, nvec, P(5,:), '-xk')
legend('t = 0', 't = 0.1 / (k_1 + k_2)', 't = 1 / (k_1 + k_2)', 't = 2 / (k_1 + k_2)', 't = 10 / (k_1 + k_2)','location','West');
%  set(l,'Interpreter','latex')
xlabel('n (number of molecules still in the first state)')
ylabel('P(t,n)');
title(strcat('Probability distributions for N=',num2str(N),' molecules undergoing reversible isomerization, at several times.'));
outsideaxes = axes('Position',[-.9 -.9 1 1],'Visible','off','YColor','w');
hold on;
set(gcf,'CurrentAxes', outsideaxes);
text(0.0001,0.0001,'CBE 504, HW1-2b, Tom Bertalan')
print(fig4,'hw1-2b.png')

K=k1/k2;
analytical_equilbrium = zeros([1,N+1]);
for n=[2:N]
    analytical_equilbrium(n) = binopdf(N+1-n,N+1,K/(1+K));
end

null_eigenvector = transpose(null(A));
%  null_eigenvector = null_eigenvector * -1; % uncomment this for matlab; comment it out for octave

fig5 = figure(5);
ax5 = axes();
hold on;
numerical = scatter(ax5, nvec, P(6,:),[],[0 0 0]);
scaled = scatter(ax5, nvec, null_eigenvector*max(analytical_equilbrium)/max(null_eigenvector), 20, [1 0 0],'+');
analytical = plot(ax5, nvec, analytical_equilbrium, 'r');
title(strcat('long-time (equilibrium) distribution for N=',num2str(N),' molecules undergoing reversible isomerization.'));
%  legend([numerical, scaled, analytical],'numerical','scaled null eigenvector','analytical','location','West') % uncomment this for matlab; comment it out for octave
xlabel('n (number of molecules still in the first state)')
ylabel(ax5, 'P(t,n)');
outsideaxes5 = axes('Position',[-.9 -.9 1 1],'Visible','off','YColor','w');
hold on;
set(gcf,'CurrentAxes', outsideaxes5);
text(0.0001,0.0001,'CBE 504, HW1-2c, Tom Bertalan')
print(fig5,'hw1-2c.png');

irreversible.m

function A=irreversible(N)
    % Returns a square [N+1, N+1] matrix for an irreversible chemical mater equation, with rate konstant k=1 s^-1
    A = zeros([N+1, N+1]);
    for i=[1:N+1]
        A(i,i) = -1*(i-1);
        if i<=N;
            A(i,i+1) = i;
        end
    end

mastereqn.m

function dpdt = mastereqn(P,A)
    dpdt = A * P;

mx_exp_solve.m

function x=mx_exp_solve(A,t,x0)
    x = expm(A*t)*x0;

ode45.m

function [tout,xout] = ode45(FUN,tspan,x0,pair,ode_fcn_format,tol,trace,count)

% Copyright (C) 2001, 2000 Marc Compere
% This file is intended for use with Octave.
% ode45.m 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.
%
% ode45.m 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 at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------
%
% ode45 (v1.11) integrates a system of ordinary differential equations using
% 4th & 5th order embedded formulas from Dormand & Prince or Fehlberg.
%
% The Fehlberg 4(5) pair is established and works well, however, the
% Dormand-Prince 4(5) pair minimizes the local truncation error in the
% 5th-order estimate which is what is used to step forward (local extrapolation.)
% Generally it produces more accurate results and costs roughly the same
% computationally.  The Dormand-Prince pair is the default.
%
% This is a 4th-order accurate integrator therefore the local error normally
% expected would be O(h^5).  However, because this particular implementation
% uses the 5th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 5th-order estimate should yield errors on the order of O(h^6).
%
% The order of the RK method is the order of the local *truncation* error, d.
% The local truncation error is defined as the principle error term in the
% portion of the Taylor series expansion that gets dropped.  This portion of
% the Taylor series exapansion is within the group of terms that gets multipled
% by h in the solution definition of the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to
% within O(h^(p+1)).  The difference between two different-order solutions is
% the definition of the "local error," l.  This makes the local error, l, as
% large as the error in the lower order method, which is the truncation
% error, d, times h, resulting in O(h^(p+1)).
% Summary:   For an RK method of order-p,
%            - the local truncation error is O(h^p)
%            - the local error (used for stepsize adjustment) is O(h^(p+1))
%
% This requires 6 function evaluations per integration step.
%
% Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableu in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
%
% Usage:
%         [tout, xout] = ode45(FUN, tspan, x0, pair, ode_fcn_format, tol, trace, count)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% pair  - flag specifying which integrator coefficients to use:
%            0 --> use Dormand-Prince 4(5) pair (default)
%            1 --> use Fehlberg pair 4(5) pair
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode45
%
% OUTPUT:
% tout  - Returned integration time points (column-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% CompereM@asme.org
% created : 06 October 1999
% modified: 17 January 2001

if nargin < 8, count = 0; end
if nargin < 7, trace = 0; end
if nargin < 6, tol = 1.e-6; end
if nargin < 5, ode_fcn_format = 0; end
if nargin < 4, pair = 0; end

pow = 1/6; % see p.91 in the Ascher & Petzold reference for more infomation.

if (pair==0),
   % Use the Dormand-Prince 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/5;
    a_(3,1)=3/40; a_(3,2)=9/40;
    a_(4,1)=44/45; a_(4,2)=-56/15; a_(4,3)=32/9;
    a_(5,1)=19372/6561; a_(5,2)=-25360/2187; a_(5,3)=64448/6561; a_(5,4)=-212/729;
    a_(6,1)=9017/3168; a_(6,2)=-355/33; a_(6,3)=46732/5247; a_(6,4)=49/176; a_(6,5)=-5103/18656;
    a_(7,1)=35/384; a_(7,2)=0; a_(7,3)=500/1113; a_(7,4)=125/192; a_(7,5)=-2187/6784; a_(7,6)=11/84;
    % 4th order b-coefficients
    b4_(1,1)=5179/57600; b4_(2,1)=0; b4_(3,1)=7571/16695; b4_(4,1)=393/640; b4_(5,1)=-92097/339200; b4_(6,1)=187/2100; b4_(7,1)=1/40;
    % 5th order b-coefficients
    b5_(1,1)=35/384; b5_(2,1)=0; b5_(3,1)=500/1113; b5_(4,1)=125/192; b5_(5,1)=-2187/6784; b5_(6,1)=11/84; b5_(7,1)=0;
    for i=1:7,
     c_(i)=sum(a_(i,:));
    end
else, % pair==1 so use the Fehlberg 4(5) coefficients:
    a_(1,1)=0;
    a_(2,1)=1/4;
    a_(3,1)=3/32; a_(3,2)=9/32;
    a_(4,1)=1932/2197; a_(4,2)=-7200/2197; a_(4,3)=7296/2197;
    a_(5,1)=439/216; a_(5,2)=-8; a_(5,3)=3680/513; a_(5,4)=-845/4104;
    a_(6,1)=-8/27; a_(6,2)=2; a_(6,3)=-3544/2565; a_(6,4)=1859/4104; a_(6,5)=-11/40;
    % 4th order b-coefficients (guaranteed to be a column vector)
    b4_(1,1)=25/216; b4_(2,1)=0; b4_(3,1)=1408/2565; b4_(4,1)=2197/4104; b4_(5,1)=-1/5;
    % 5th order b-coefficients (also guaranteed to be a column vector)
    b5_(1,1)=16/135; b5_(2,1)=0; b5_(3,1)=6656/12825; b5_(4,1)=28561/56430; b5_(5,1)=-9/50; b5_(6,1)=2/55;
    for i=1:6,
     c_(i)=sum(a_(i,:));
    end
end % if (pair==0)


% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/1e9;
h = (tfinal - t)/100; % initial guess at a step size
x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0; end
end % if count

if trace
 clc, t, h, x
end

if (pair==0),

   k_ = x*zeros(1,7);  % k_ needs to be initialized as an Nx7 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)

   % Compute the first stage prior to the main loop.  This is part of the Dormand-Prince pair caveat.
   % Normally, during the main loop the last stage for x(k) is the first stage for computing x(k+1).
   % So, the very first integration step requires 7 function evaluations, then each subsequent step
   % 6 function evaluations because the first stage is simply assigned from the last step's last stage.
   % note: you can see this by the last element in c_ is 1.0, thus t+c_(7)*h = t+h, ergo, the next step.
   if (ode_fcn_format==0), % (default)
      k_(:,1)=feval(FUN,t,x); % first stage
   else, % ode_fcn_format==1
      k_(:,1)=feval(FUN,x,t);
   end % if (ode_fcn_format==1)

   % increment rhs_counter
   if (count==1), rhs_counter = rhs_counter + 1; end

   % The main loop using Dormand-Prince 4(5) pair
   while (t < tfinal) & (h >= hmin),
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nxs, a_ is 7x6, which is s by (s-1),
      %        s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages)
      if (ode_fcn_format==0), % (default)
         for j = 1:6,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:6,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_*b4_); % k_ is Nxs (or Nx7) and b4_ is a 7x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is Nxs (or Nx7) and b5_ is a 7x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta<=tau),
         t = t + h;
         x = x5;    % <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end % if (delta<=tau)

      if trace
         home, t, h, x
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

      % Assign the last stage for x(k) as the first stage for computing x(k+1).
      % This is part of the Dormand-Prince pair caveat.
      % k_(:,7) has already been computed, so use it instead of recomputing it
      % again as k_(:,1) during the next step.
      k_(:,1)=k_(:,7);

   end % while (t < tfinal) & (h >= hmin)

else, % pair==1 --> use RK-Fehlberg 4(5) pair

   k_ = x*zeros(1,6);  % k_ needs to be initialized as an Nx6 matrix where N=number of rows in x
                       % (just for speed so octave doesn't need to allocate more memory at each stage value)
 
   % The main loop using Dormand-Prince 4(5) pair
   while (t < tfinal) & (h >= hmin),
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
      % notes: k_ needs to end up as an Nx6, a_ is 6x5, which is s by (s-1),  (RK-Fehlberg has s=6 stages)
      %        s is the number of intermediate RK stages on [t (t+h)]
      if (ode_fcn_format==0), % (default)
     k_(:,1)=feval(FUN,t,x); % first stage
      else, % ode_fcn_format==1
     k_(:,1)=feval(FUN,x,t);
      end % if (ode_fcn_format==1)

      if (ode_fcn_format==0), % (default)
         for j = 1:5,
            k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)');
         end
      else, % ode_fcn_format==1
         for j = 1:5,
            k_(:,j+1) = feval(FUN, x+h*k_(:,1:j)*a_(j+1,1:j)', t+c_(j+1)*h);
         end
      end % if (ode_fcn_format==1)

      % increment rhs_counter
      if (count==1), rhs_counter = rhs_counter + 6; end

      % compute the 4th order estimate
      x4=x + h* (k_(:,1:5)*b4_); % k_(:,1:5) is an Nx5 and b4_ is a 5x1

      % compute the 5th order estimate
      x5=x + h*(k_*b5_); % k_ is the same as k_(:,1:6) and is an Nx6 and b5_ is a 6x1

      % estimate the local truncation error
      gamma1 = x5 - x4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(x,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if (delta<=tau),
         t = t + h;
         x = x5;    % <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end % if (delta<=tau)

      if trace
         home, t, h, x
      end % if trace

      % Update the step size
      if (delta==0.0),
       delta = 1e-16;
      end % if (delta==0.0)
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end % while (t < tfinal) & (h >= hmin)


end % if (pair==0),

if (t < tfinal)
   disp('Step size grew too small.')
   t, h, x
end

reversible.m

function A=reversible(N,kvec)
    % Returns a square [N+1, N+1] matrix for a reversible chemical master equation, with rate konstants given in the vector kvec.
    A=zeros([N+1,N+1]);
    k1=kvec(1);
    k2=kvec(2);
    %  first row:
    n=0;
    A(1,2)=k1;
    A(1,1)=-1*k2*N;
    %  Last row:
    n=N;
    A(N+1,n+1)=-1*k1*N;
    A(N+1,N)=k2;
    %  Other rows:
    for i=[2:N]
        n=i-1;
        A(i,i)   = -1*(k1*n + k2*(N-n));
        A(i,i-1) = k2*(N-n+1);
        A(i,i+1) = k1*(n+1); 
    end