% % NRM_LV_timetrack.m Run Gillespie's algorithm on a Lotka Volterra model. % % usage: GILL_LV_timetrack(tend) % % tend = max time for the program to run = 20 % % example: GILL_LV_timetrack(20) % function [X] = GILL_LV_timetrack(tend) %maximum number of steps to run before quitting. maxi = tend*10^5; %rate constants. k = [10, 0.01, 0.01, 0.01, 10]; %reaction vectors (columns); zeta= [1 0 -1 -1 0; 0 1 0 0 -1]; %initialize. Assume start with 300 animals each. initial = [300,300]'; %T is time. T = zeros([maxi,1]); %state of the system. x = zeros([2,maxi]); %initialize the state of the system x(:,1) = initial; %holder for intensities. lambda = zeros([5,1]); %main loop for i =1:maxi-1 %set intensity functions. lambda(1) = k(1)*x(1,i); lambda(2) = k(2)*x(1,i)*x(2,i); lambda(3) = k(3)*x(1,i)*x(2,i); lambda(4) = k(4)*x(1,i); lambda(5) = k(5)*x(2,i); %find the sum of the intensity functions lambda0 = lambda(1) + lambda(2) + lambda(3) + lambda(4) + lambda(5); %generate a RV for finding when the next reaction takes place. r1 = rand; %find the time of the occurrence. t = log(1/r1)/lambda0; %if the time takes us past the end time, we should quit. if T(i) + t > tend T(i+1) = tend; x(:,i+1) = x(:,i); break end %find which reaction occured. r2 = rand; if r2 < lambda(1)/lambda0 loc = 1; elseif r2 < (lambda(1) + lambda(2))/lambda0 loc = 2; elseif r2 < (lambda(1) + lambda(2) + lambda(3))/lambda0 loc = 3; elseif r2 < (lambda(1) + lambda(2) + lambda(3) + lambda(4))/lambda0 loc = 4; else loc = 5; end %update the state of the system and the time. x(:,i+1) = x(:,i) + zeta(:,loc); T(i+1) = T(i) + t; end count = i; X = x(count+1); %plot the timecourse of the state of the system. figure plot(T(1:count),x(1,(1:count)),T(1:count),x(2,(1:count))) legend('Rabbits','Foxes') % title('Full version') %Now for the exact solution % % Num = 1000000; % Delta = tend/Num; % % ex = zeros([2,Num+1]); % ex(:,1) = initial; % T = zeros([Num+1,1]); % for i = 1:Num % lambda(1) = k(1)*ex(1,i); % lambda(2) = k(2)*ex(1,i)*ex(2,i); % lambda(3) = k(3)*ex(1,i)*ex(2,i); % lambda(4) = k(4)*ex(1,i); % lambda(5) = k(5)*ex(2,i); % % ex(:,i+1) = ex(:,i) + zeta*lambda*Delta; % T(i+1) = T(i) + Delta; % end % % figure % plot(T,ex(1,:),T,ex(2,:)) % legend('Rabbits','Foxes') end %The program