Gillespie算法
问题描述：它是一个随机模拟算法，在这里有四个反应，三个反应物中，在一个小区间[t,t+tau]内每次只能发生一个反应,通过产生两个随机数，再根据公式来更新时间个状态疑问：程序不能运行到规定时间就暂停了，但是标准答案是可以运行到规定的结束时间，通过观察运行出来的数据发现是因为的a00=0，导致tau=inf（一般情况下，tau是很小的数），导致程序直接结束
希望达到的效果：一直更新tau，直至到规定的截止时间
代码：
clear;
clc;
c1=0.1;c2=0.02;c3=0.5;c4=0.04 ; %rate constant
v=[-1 -2 2 0;0 1 -1 -1;0 0 0 1]; %update matrix
count=1;init=[100;0;0]; s(:,count)=init;t(count)=0;%initial
length=400; %length of time for computation
xold=s(:,count);told=t(count);
while told<length
r1=rand;r2=rand;
a(1)=c1*xold(1);a2=c2*xold(1)*(xold(1)-1)/2;a(3)=c3*xold(2);a(4)=c4*xold(2);
a00=sum(a);
tau=(1/a00)*log(1/r1); %find stepsize tau
sss=0;kkk=0;judge='fff';standard=r2*a00; %find reaction
while strcmp(judge,'fff')
kkk=kkk+1;sss=sss+a(kkk);
if sss>=standard
index=kkk;judge='ttt';
break;
end
end %reaction #:index
told=told+tau;xold=xold+v(:,index); %update
count=count+1;t(count)=told;s(:,count)=xold; %save state to t&s
end %for the do while for one simulation
plot(t,s(1,:),'r');hold on
plot(t,s(2,:),'g');hold on
plot(t,s(3,:),'y');
[此贴子已经被作者于2022-5-18 10:21编辑过]
