% A program for exit calculations for MBI-NIMBioS-CAMBAM PROGRAM % Euler's method is used for solving the SDEs % nsamp sample paths are calculated % The proportion exiting up to time t are calculated % Exits occur when either population size is less than 1 % y1 and y2 are the two populations (different species) % b1,d1,b2,d2 are the per capita birth and death rates % y10 y20 are the initial population sizes % Problem-dependent statements are marked with a %*** (9 statements) % Accuracy generally increases as h decreases % icase=1 corresponds to the deterministic problem and so nsamp=1 clf clear diary outputsde.m for icase=1:2 clear tt clear yp1 clear yp2 nsamp=1000; %*** tmax=100; %*** nt=1000; %*** y10=15; %*** y20=15; %*** if(icase==1) nsamp=1; end h=tmax/nt; hs=sqrt(h); randn('state',2); %initiates the random number generator te1=zeros(nsamp,1); te2=zeros(nsamp,1); te3=zeros(nsamp,1); jj1=0; jj2=0; jj3=0; for jj=1:nsamp y1=y10; y2=y20; yp1(1)=y1; yp2(1)=y2; r=randn(nt+1,2); nchk1=0; nchk2=0; n=0; t=0; chk=0; tt(1)=0; while (chk==0) n=n+1; t=t+h; if(jj==nsamp) tt(n+1)=t; end b1=.84; %*** d1=.40+.01*y1+.022*y2; %*** b2=.90; %*** d2=.75+.0067*y2+.005*y1; %*** f1=b1*y1-d1*y1; f2=b2*y2-d2*y2; g1=sqrt(b1*y1+d1*y1); g2=sqrt(b2*y2+d2*y2); if(icase==1) g1=0; end if(icase==1) g2=0; end y1=y1+h*f1+hs*g1*r(n,1); y2=y2+h*f2+hs*g2*r(n,2); if(jj==nsamp) yp1(n+1)=y1; end if(jj==nsamp) yp2(n+1)=y2; end % This is Euler's approximation to the SDE if (y1 < 1) chk=1; jj1=jj1+1; te1(jj1)=t; end if (y2 < 1) chk=1; jj2=jj2+1; te2(jj2)=t; end if (t > tmax) chk=1; jj3=jj3+1; te3(jj3)=t; chk=1; end end % end of while (chk==0) loop end % end of for jj=1:nsamp loop tp=0; tp1=0; tp2=0; tp3=0; if(jj1 ~= 0) tp1=sum(te1)/jj1; end if(jj2 ~= 0) tp2=sum(te2)/jj2; end if(jj3 ~= 0) tp3=sum(te3)/jj3; end if(jj1+jj2 ~= 0) tp=(sum(te1)+sum(te2))/(jj1+jj2); end p1=jj1/nsamp; p2=jj2/nsamp; p3=jj3/nsamp; % tp1 and tp2 are mean exit times for populations 1 and 2 % tp3 is mean time for sample paths not exiting % tp is mean exit time for paths that exit % p1, p2 are proportions exiting for populations 1 and 2 % p3 is proportion not exiting in time tmax disp(' ') if(icase==1) disp(' Deterministic Calculational Results'); end if(icase==2) disp(' Stochastic Calculation Results'); end disp(' icase nsamp h tmax') disp((sprintf(' %12.0f %12.0f %12.5f %12.2f',icase,nsamp,h,tmax))); disp(' tp1 p1') disp((sprintf(' %12.6f %12.6f', tp1, p1))); disp(' tp2 p2') disp((sprintf(' %12.6f %12.6f', tp2, p2))); disp(' tp3 p3') disp((sprintf(' %12.6f %12.6f', tp3, p3))); disp(' tp p1+p2') disp((sprintf(' %12.6f %12.6f', tp, p1+p2))); subplot(2,2,2*icase-1) set(gca,'fontsize',18,'linewidth',1.5); plot(yp1,yp2,'k-') xlabel('Pop. y1') ylabel('Pop. y2') if(icase==1) title('Deterministic'); end %if(icase==2) title('Stochastic'); end hold on subplot(2,2,2*icase) set(gca,'fontsize',18,'linewidth',1.5); plot(tt,yp1,'r-',tt,yp2,'k-') xlabel('Time t') ylabel('Pops. 1 and 2') if(icase==1) title('Deterministic'); end %if(icase==2) title('Stochastic'); end hold on end % end of for icase=1:2 loop hold off diary off