% % NRM_LV_timetrack.m Run NRM on a Lotka Volterra model. % % usage: NRM_LV_timetrack(tend) % % tend = max time for the program to run = 20 % % example: NRM_LV_timetrack(20) % function [] = NRM_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 real 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]); %Set Poisson times and next Poisson Jumps. Also set t's. Tk = zeros([5,1]); Pk = zeros([5,1]); t = zeros([5,1]); %set the first exponential holding times of the Poisson processes. for i = 1:5 r = rand; Pk(i) = log(1/r); end %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); %set t's, the times until each of the reactions would go if left alone. % This is equivalent to finding how long each "exponential alarm clock" % would require. for ell = 1:5 if lambda(ell) > 0 t(ell) = (Pk(ell) - Tk(ell))/lambda(ell); if t(ell)<0 t end else t(ell) = inf; end end %loc gives the index of the minimum and t(loc) is the time until the %first alarm clock would go off. loc = 1; for ind = 2:5 if t(loc)>t(ind) loc = ind; end end %if the time takes us past the end time, we should quit. if T(i) + t(loc)> tend T(i+1) = tend; x(:,i+1) = x(:,i); break end %update the state of the system and the time. x(:,i+1) = x(:,i) + zeta(:,loc); T(i+1) = T(i) + t(loc); %update the internal times of the Poisson processes. for m = 1:5 Tk(m) = Tk(m) + lambda(m)*t(loc); end %We used one exponential random variable (the one from the Poisson %process that fired). We need to set the next firing time of that %Poisson process. r = rand; Pk(loc) = Pk(loc) + log(1/r); end count = i; %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