8.2 Programs
147
v(i+1,j) = v(i-1,j) + (-dGamma*v(i,j)- Dd1U(i,j) -0.5*Dlambda*Dd3U(i,j) +. . .
F(i) + dGamma*F1(i) + Fext )*2.*dt;
% to add in case you want to simulate an external force Fext = eta(i,j)*FDN, tau = ;
%========================================================
% NOISE
%========================================================
a1 = abs(randn(1));
b1 = abs(randn(1));
a2 = abs(randn(1));
b2 = abs(randn(1));
a3 = abs(randn(1));
b3 = abs(randn(1));
h1(i,j) = sqrt((-2.*(DDeff(i,j)*D1)/(t1*tR))*(1-E1 2 )*log(a1))*cos(2*pi*b1);
h2(i,j) = sqrt((-2.*(DDeff(i,j)*D2)/(t2*tR))*(1-E2 2 )*log(a2))*cos(2*pi*b2);
h3(i,j) = sqrt((-2.*(DDeff(i,j)*D3)/(t3*tR))*(1-E3 2 )*log(a3))*cos(2*pi*b3);
eps1(i+1,j)=eps1(i,j)*E1 + h1(i,j);
eps2(i+1,j)=eps2(i,j)*E2 + h2(i,j);
eps3(i+1,j)=eps3(i,j)*E3 + h3(i,j);
eps(i+1,j)=(eps1(i+1,j)+eps2(i+1,j)+eps3(i+1,j));
%=========================================================
end
end
MEANXR=(sum(x,2))/NR;
MEANVR=(sum(v,2))/NR;
MEANEPS=(sum(eps,2))/NR;
MEANDeff=(sum(DDeff,2))/NR;
t=(1:N)*dt;
t=t’;
save t.dat t -ascii;
figure(1); plot(t,MEANXR); grid;
figure(2); plot(t,MEANVR); grid;
figure(3); plot(MEANXR,MEANVR); grid;
save XCg1E-1-thetapi.DAT MEANXR -ascii;
save VCg1E-4-thetapi.DAT MEANVR -ascii;
%=========================================================