% A two-patch SIS program for MBI-NIMBioS-CAMBAM PROGRAM % Euler's method is used for solving the SDEs % nsamp sample paths are calculated % s1,s2,i1,i2 are population sizes in two patches % gam1,gam2,di12,di21,ds12,ds21,bet1,bet2 are constants % s10,s20,i10,i20 are the initial population sizes % Problem-dependent statements are marked with a %*** (15 statements) % Accuracy generally increases as h decreases % icase=1 corresponds to the deterministic problem and so nsamp=1 clf clear %diary outputgroupsde.m for icase=1:2 clear tt clear s1 clear s2 clear i1 clear i2 nsamp=1000; %*** tmax=25; %*** nt=500; %*** s10=25; %*** s20=20; %*** i10=30; %*** i20=25; %*** gam1=.05; %*** gam2=.4; %*** di12=.3; %*** di21=.05; %*** ds12=0; %*** ds21=.25; %*** bet1=.5; %*** bet2=.01; %*** if(icase==1) nsamp=1; end h=tmax/nt; hs=sqrt(h); randn('state',3); %initiates the random number generator ss1=0; ss12=0; ss2=0; ss22=0; si1=0; si12=0; si2=0; si22=0; for jj=1:nsamp s1=s10; s2=s20; i1=i10; i2=i20; s1p(1)=s1; s2p(1)=s2; i1p(1)=i1; i2p(1)=i2; r=randn(nt+1,8); 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 nt1=s1+i1; nt2=s2+i2; fs1=-s1*i1*bet1/nt1+gam1*i1+ds21*s2-ds12*s1; fi1=s1*i1*bet1/nt1-gam1*i1-di12*i1+di21*i2; fs2=-s2*i2*bet2/nt2+gam2*i2-ds21*s2+ds12*s1; fi2=s2*i2*bet2/nt2-gam2*i2-di21*i2+di12*i1; hlp1=-sqrt(s1*i1*bet1/nt1)*r(n,1)+sqrt(gam1*i1)*r(n,2); hlp2=sqrt(ds21*s2)*r(n,3)-sqrt(ds12*s1)*r(n,4); hlp3=sqrt(di21*i2)*r(n,5)-sqrt(di12*i1)*r(n,6); hlp4=-sqrt(s2*i2*bet2/nt2)*r(n,7)+sqrt(gam2*i2)*r(n,8); gs1= hlp1+hlp2; gi1= -hlp1+hlp3; gs2= hlp4-hlp2; gi2= -hlp4-hlp3; if(icase==1) gs1=0; end if(icase==1) gs2=0; end if(icase==1) gi1=0; end if(icase==1) gi2=0; end s1=s1+h*fs1+hs*gs1; i1=i1+h*fi1+hs*gi1; s2=s2+h*fs2+hs*gs2; i2=i2+h*fi2+hs*gi2; % This is Euler's approximation to the SDE if(jj==nsamp) s1p(n+1)=s1; end if(jj==nsamp) i1p(n+1)=i1; end if(jj==nsamp) s2p(n+1)=s2; end if(jj==nsamp) i2p(n+1)=i2; end %if (s1 < 0) % i1=i1+s1; % s1=0; %end %if (i1 < 0) % s1=s1+i1; % i1=0; %end %if (s2 < 0) % i2=i2+s2; % s2=0; %end %if (i2 < 0) % s2=s2+i2; % i2=0; %end if (t > tmax) ss1=s1+ss1; ss12=s1*s1+ss12; ss2=s2+ss2; ss22=s2*s2+ss22; si1=i1+si1; si12=i1*i1+si12; si2=i2+si2; si22=i2*i2+si22; chk=1; end end % end of while (chk==0) loop end % end of for jj=1:nsamp loop ss1=ss1/nsamp; ss12=ss12/nsamp-ss1*ss1; ss2=ss2/nsamp; ss22=ss22/nsamp-ss2*ss2; si1=si1/nsamp; si12=si12/nsamp-si1*si1; si2=si2/nsamp; si22=si22/nsamp-si2*si2; ntot=ss1+ss2+si1+si2; disp(' ') if(icase==1) disp(' ODE Calculational Results'); end if(icase==2) disp(' SDE Calculational Results (Means and Standard Deviations)'); end disp(' icase nsamp h tmax') disp((sprintf(' %12.0f %12.0f %12.5f %12.2f',icase,nsamp,h,tmax))); disp(' S1(mean) I1(mean) S2(mean) I2(mean)') disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', ss1, si1, ss2, si2))); disp(' S1(std) I1(std) S2(std) I2(std)') disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', sqrt(ss12), sqrt(si12), sqrt(ss22), sqrt(si22)))); disp(' ntot ') disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', ntot))); subplot(2,2,2*icase-1) plot(tt,s1p,'k-',tt,i1p,'r--','linewidth',2) xlabel(' Time t') ylabel('S_1 and I_1') legend('S_1', 'I_1'); axis([0,tmax,0,60]) if(icase==1) title('ODE'); end if (icase==2) title('SDE'); end hold off subplot(2,2,2*icase) plot(tt,s2p,'k-',tt,i2p,'r--','linewidth',2) xlabel('Time t') ylabel('S_2 and I_2') legend('S_2', 'I_2') axis([0,tmax,0,60]) if(icase==1) title('ODE'); end if(icase==2) title('SDE'); end hold off end % end of for icase=1:2 loop hold off %diary off