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:27编辑过]