Figure 5.15:

Cumulative distribution for 2 A \mathrel {\mkern 4mu}\mathrel {\smash {\mathchar "392D}}\mathrel {\mkern -2.5mu}-> B at t=1 with n_0=500, \Omega =500. Discrete master equation (steps) versus omega expansion (smooth).

Code for Figure 5.15

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
% solve the pde for the random walk and compare the
% discrete chemical master equation
% jbr, 12/25/2011
%
% passed parameters to functions, jbr, 6/27/2018
%

omega = 500;
k = 1;
c0 = 1;
n0 = c0*omega;
b = 3;

ncol = 25;
[x, A, B, q] = colloc (ncol-1, 'left');
x = b*x;
A = A/b;
B = B/(b*b);
q = q*b;

p = struct();
p.k = k; p.c0 = c0; p.x = x; p.A = A; p.B = B;


tstop = 1;
nts = 25;
tsteps = linspace(0, tstop, nts)';
var = 1e-2;
% y0 = ones(ncol,1);
% y0(1) = 1/2;
y0 = normcdf(x, 0, 5*var);
ydot0 = zeros(ncol,1);
resid = continuous(0, y0, ydot0, p);
ydot0 = -resid;
ydot0(1) = 0;

opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tout, y] = ode15i (@(t, F, Fdot) continuous(t, F, Fdot, p), tsteps, y0, ydot0, opts);

cstop = 1./(1/c0+k*tstop);
xshift = cstop + [-flipud(x); x]/sqrt(omega);
F = [1-fliplr(y(end,:)), y(end,:)]';

%compare master equation

n = (0:n0)';
d1 = -k*n.*(n-1)/(2*omega);
d2 = k*(n+2).*(n+1)/(2*omega);
d2(end-1:end) = [];
%create A matrix with d1 on diagonal and d2 on 2nd super-diagonal
A = diag(d1) + diag(d2,2);
% initialize with n0 A molecules
P0 = zeros(length(n),1);
P0(end) = 1;
% solve master equation at time t
t = 1;
P = expm(A*t)*P0;
Fmas = cumsum(P);

nscal = n/omega;

figure()
plot(nscal, Fmas, xshift, F)
axis ([0.4,0.6,0,1])

figure()
plot(nscal, Fmas, xshift, F)

results1 = [nscal, Fmas];
results2 = [xshift, F];

save F2nd.dat results1 results2

continuous.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
function resid = continuous(t, F, Fdot, p)
  c = 1/(1/p.c0+p.k*t);
  v = -2*p.k*c*p.x;
  D=p.k*c*c;
  Fx = p.A*F;
  dFdt = -v.*Fx + D*p.B*F;
  % pde
  resid = Fdot - dFdt;
  % F = 1/2 at x =0; F=1 at large x
  resid(1) = F(1) - 0.5;
  resid(end) = F(end) -1;