Figure 5.14:

Simulation of 2 A \mathrel {\mkern 4mu}\mathrel {\smash {\mathchar "392D}}\mathrel {\mkern -2.5mu}-> B for n_0=500, \Omega =500. Top: discrete simulation; bottom: SDE simulation.

Code for Figure 5.14

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
% Run KMC samples for  2A -> B
% Compare SDE approximation
%
% jbr, 12/20/2011
%

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

nsim = 10*omega;
time = zeros (nsim+1, 1);
time(1) = 0;
n = zeros(nsim+1,1);
n(1) = n0;
rand('seed', -10);
for i=1:nsim
  r = k*n(i)*(n(i)-1)/(2*omega);
  if (r == 0)
    n = n(1:i); time = time(1:i);
    break
  end
  % choose time
  tau = -log(rand)/r;
  time(i+1)=time(i)+tau;
  % fire reaction
  n(i+1) = n(i) - 2;
end

% deterministic solution
nts = 500;
tinit = 1e-2;
tfin = 1e2;
tc = logspace(log10(tinit), log10(tfin), nts)';
c = 1./(1/c0 + k*tc); c = c(:);

% sde solution at time points of deterministic solution
xi = zeros(nts,1);
% rand('seed',0);
for i = 1:nts-1
  Delta = tc(i+1)-tc(i);
  xi(i+1) = xi(i) -2*k*c(i)*xi(i)*Delta + sqrt(2*k*c(i)*c(i)*Delta)*randn;
end
csde = c + xi/sqrt(omega);

figure()
[ts, ns] = stairs(time, n);
cs = ns/omega;
semilogx(ts, cs, tc, c)
axis ([tinit, tfin, 0, 1])

figure()
semilogx(tc, csde, tc, c)
axis ([tinit, tfin, 0, 1])

data = [ts(1:end-1), cs(1:end-1), tc, c, csde];
save "sam2nd.dat" data