42
3 Biased Brownian Motion
3.7 Programs
3.7.1 Program 3.1, Euler Equation, Matlab Code
%Solve the Euler Equation
N=20000; NR=1000; dt=0.01; D=0.0005; Vo=1.; F=0.; x00=0.;
raizdt=sqrt(2*D*dt);
dw=raizdt*randn(N,NR);
x=zeros(N-1,NR);
V=zeros(N-1,NR);
dV=zeros(N-1,NR);
for j=1:NR
x(1,j)=x00;
for i=1:N-1
xij=x(i,j);
%V(i,j)=(Vo/(2*pi))*(sin(2*pi*xij/L)- 0.25*sin(4*pi*xij/L));
%dV(i,j)=(Vo/(L))*(cos(2*pi*xij/L)- 0.5*cos(4*pi*xij/L));
%dV is the derivative of V, “i” is the time and “j” is the stochastic realization
dV(i,j)=(cos(2.*pi*xij) - 0.5*cos(4.*pi*xij)); %dimensionless
x(i+1,j)=xij - Vo*dV(i,j)*dt + D*F*dt + dw(i,j);%dimensionless
end
end
MEANXR=(sum(x,2))/NR;%Mean over the stochastic realizations.
t=((1:N)*dt)’;
figure(1); plot(t,MEANXR); grid;
save t.dat t -ascii; save MEANXR.DAT MEANXR -ascii;
3.7.2 Program 3.2, F-P Equation, Matlab Code
%Solution of the Dimensionless Fokker-Planck Equation by Crank-Nicholson
scheme
%V(x)=(1/2pi)[sin(2pix) - 0.25sin(4pix)], Potential slightly tilted to the right
%Intensity s0, Central Position x0, Standard deviation sigma,
%Diffusion coefficient D
%V (x) = d1, V "(x) = d2,
x0=0.5; N=100; tmax=100; dt=0.01; dx=0.01; D=0.0005; F=6.; Vo=1; sigma=0.1
k=0.5*D*dt/dx 2 ;
s2=1/(2*(sigma 2 ));
s0=1./(sigma*sqrt(2*pi)); T=floor(tmax/dt)+1; %T is the number of temporal points
x=((1:N)*dx)’; % spatial grade, column vector
p=s0*exp(-((x − x0). 2 )*s2);%Initial Probability
p1=p;
E=zeros(N); M=zeros(N); pt=ceil(T/8);P=zeros(N,8);
for j=1:N;