62
4 The Smoluchowski Model
E1(j,j-1)=-d1;
E1(j-1,j)=-c1;
L1(j,j-1)=d1;
L1(j-1,j)=c1;
end
%STATE 2, DIAGONAL NEIGBOURGH
for j=N+2:2*N
d12=0.;
c2=ka*(1.+ka2*d12);
d2=ka*(1.-ka2*d12);
E2(j,j-1)=-d2;
E2(j-1,j)=-c2;
L2(j,j-1)=d2;
L2(j-1,j)=c2;
end
for j=N+1:2*N
KR21(j-N,j)=k21*dt;
KL12(j,j-N)=k12*dt;
end
E=E+E1+E2;
M = M + L1-K12 + L2-K21 + KR21 +KL12;
IE=inv(E);
pt=ceil(T/8);
for tp=5:5:25
for it=1:pt
f=(IE)*(M*f);
end
F(:,tp)=f; % distribution at times tp
end
figure(1); plot(x,F(:,5),x,F(:,10),x,F(:,15),x,F(:,20))
F1=F(:,5);
F2=F(:,10);
F3=F(:,15);
F4=F(:,20);
save F1.dat F1 -ascii;
save F2.dat F2 -ascii;
save F3.dat F3 -ascii;
save F4.dat F4 -ascii;
save x.dat x -ascii;
%=============================================================
%Flux1 IN STATE 1
%=============================================================
G4=zeros(2*N,1);
Flux12=zeros(2*N-1,1);
dV1=4.*(cos(2.*pi*x)-cos(4.*pi*x)+cos(6.*pi*x));