how to generate (multiple) realisations using Gillespie Algorithms

1 回表示 (過去 30 日間)
yamen alharbi
yamen alharbi 2020 年 12 月 24 日
コメント済み: yamen alharbi 2021 年 1 月 3 日
Hello, I would like to simulate this reaction
a->X
X-> X
at rates c. here is a code was written by Desmond Higham and I tried to adopted it to this reaction.
I would like to simulate different realisations . Could you please shed a light on this.
Your help is highly aprrecaited.
clf
clear all;
clc
rand('state',100)
%stoichiometric matrix
V = [1 -1];
%%%%%%%%%% Parameters and Initial Conditions %%%%%%%%%
nA = 6.023e23; % Avagadro's number
vol = 1e-15; % volume of system
X = zeros(1,1);
X(1) = 100; % molecules of substrate
c(1) = .4; c(2) =.1; c(3)=8; c(4)=.3;
t = 0;
tfinal = 100;
count = 1;
tvals(1) = 0;
Xvals(:,1) = X;
while t < tfinal
a(1) = c(1);% calculate propensity functions
a(2) = c(2)*X(1);
asum = sum(a);
j = min(find(rand<cumsum(a/asum))); % select which reaction
tau = log(1/rand)/asum; % select time
X = X + V(:,j); %update X
count = count + 1;
t = t + tau; %Update time
tvals(count) = t;
Xvals(:,count) = X;
end
%%%%%%%%%%% Plots %%%%%%%%
L = length(tvals);
tnew = zeros(1,2*(L-1));
tnew(1:2:end-1) = tvals(2:end);
tnew(2:2:end) = tvals(2:end);
tnew = [tvals(1),tnew];
%
Svals = Xvals(1,:);
ynew = zeros(1,2*L-1);
ynew(1:2:end) = Svals;
ynew(2:2:end-1) = Svals(1:end-1);
plot(tnew,ynew,'g-')
hold on
  5 件のコメント
Star Strider
Star Strider 2021 年 1 月 3 日
The only one that appears to be accessable is the SIAM paper. I’ll download it now and read it later. (The first is behind a paywall, the second offers only short descriptions and no links to other information.)
yamen alharbi
yamen alharbi 2021 年 1 月 3 日
Thank you so much
Here is the first paper by Gillespie
and I have uploaded a copy of the same paper
Many Thanks again

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeChemistry についてさらに検索

製品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by