urngame.m

%                                              
% Ehrnfest's Lottery Simulation               
% Hank Ashbaugh                               
% Tulane University                           
% Chemical and Biomolecular Engineering       
% 2009       
%

%           
% Load the variables for placement of balls in the urns (urn.dat)
% and for the number of lottery spins to be performed
%

%
% This section reads in the run and urn variables and initializes
% the lottery
%

clear
load urn.dat
load run.dat
numurns=size(urn,1);
save=10;
savecount=1;
urn1(1:numurns,1)=urn;
urn2(1:numurns,1)=urn;
urnold=urn;
numballs = sum(urn);
store=round(run/2000-.5);
if (store < 1) 
    store=1;
end
tol=1e-200;
ijk=1;

%
% Calculate the initial entropy
%

ent1(1)=numballs*log(numballs);
ent2(1)=numballs*log(numballs);
x(ijk)=ijk-1;
for i=1:numurns
    ent1(1)=ent1(1)-urn(i,1)*log(urn(i,1)+tol);
    ent2(1)=ent2(1)-urnold(i,1)*log(urnold(i,1)+tol);
end
ent1(1)=ent1(1)/numballs;
ent2(1)=ent2(1)/numballs;
nocp1(1)=urn(1,1);
nocp2(1)=urnold(1,1);

%
% This loop iterates the lottery spins
%

for i = 1:run
    
%
% Picks a marble/ball to move
%

    pickball = randi(numballs);
    number   = 0;

%
% Find the urn the picked marble is in
%

    for count = 1:numurns
        index = urn(count,1);
        if ((pickball > number) && (pickball <= (index+number)))
            here = count;
        end
        number = number + index;
    end

%
% Picks a direction for the marble to move in the multiurn lottery
%

    direct = randi(2);
    if (direct == 1)
        direct = - 1;
    else
        direct =   1;
    end
    
%
% If there are only two urns transfer marble to other urn
%

    if ((numurns == 2) && (here == 1)) 
        direct =  1;
    end
    if ((numurns == 2) && (here == 2)) 
        direct = -1;
    end

%
% Move the marble and solve the master equation for spin s to s+1
%

    now = here + direct;
    if ((now > 0) && (now <= numurns)) 
        urn(here,1) = urn(here,1) - 1;
        urn(now,1)  = urn(now,1)  + 1;
    end
    if (numurns == 2)
        urnnew(1,1) = urnold(1,1)+(urnold(2,1)-urnold(1,1))/numballs;
        urnnew(numurns,1) = urnold(numurns,1)+(urnold(numurns-1,1)-urnold(numurns,1))/numballs;
    else
        for count = 2:(numurns-1)
            urnnew(count,1) = urnold(count,1)+(urnold(count-1,1)+urnold(count+1,1)-2*urnold(count,1))/(2*numballs);
        end
        urnnew(1,1) = urnold(1,1)+(urnold(2,1)-urnold(1,1))/(2*numballs);
        urnnew(numurns,1) = urnold(numurns,1)+(urnold(numurns-1,1)-urnold(numurns,1))/(2*numballs);
    end
    urnold = urnnew;
    j=round(i/store-.5);

%
% Calculate the lottery entropy
%

    if ((i-j*store) == 0)
        ijk=ijk+1;
        ent1(ijk)=numballs*log(numballs);
        ent2(ijk)=numballs*log(numballs);
        x(ijk)=i;
        for k=1:numurns
         ent1(ijk)=ent1(ijk)-urn(k,1)*log(urn(k,1)+tol);
         ent2(ijk)=ent2(ijk)-urnold(k,1)*log(urnold(k,1)+tol);
        end
        ent1(ijk)=ent1(ijk)/numballs;
        ent2(ijk)=ent2(ijk)/numballs;
%
% For the two urn game periodically save urn occupations 
% for plotting at the end of the run
%
        if (numurns == 2)
            nocp1(ijk)=urn(1,1);
            nocp2(ijk)=urnold(1,1);
        end
    end
    
%
% For the multi-urn game save urn occupations on the 10^i (i = 1,2,...)
% spin for plotting at the end of the run
%

    j=round(i/save-.5);
    if ((i-j*save) == 0)
        savecount=savecount+1;
        save=save*10;
        urn1(1:numurns,savecount)=urn;
        urn2(1:numurns,savecount)=urnold;
    end
end

%
% Plot urn occupations for the two and multi-urn lotteries
%

if (numurns == 2)
    uu(1:ijk,1)=nocp1;
    uu(1:ijk,2)=nocp2;
    uu(1:ijk,3)=numballs/2;
    subplot(2,1,1)
    plot(x,uu)
    axis([1,x(ijk),0,numballs])
    title('random and master equation lotteries')
    xlabel('spin')
    ylabel('n1 : urn 1 occupation')
else
    subplot(2,2,1)
    plot(urn1)
    axis([1,numurns,0,inf])
    title('random lottery ')
    xlabel('urn')
    ylabel('ni : urn i occupation')
    subplot(2,2,2)
    plot(urn2)
    axis([1,numurns,0,inf])
    title('master equation lottery')
    xlabel('urn')
    ylabel('<ni> : urn i mean occupation')
end

%
% Plot the entropy for the two and multi-urn lotteries
%

entropy(1:ijk,1)=ent1(1:ijk);
entropy(1:ijk,2)=ent2(1:ijk);
entropy(1:ijk,3)=log(numurns);
subplot(2,1,2)
plot(x,entropy)
axis([1,x(ijk),0,1.1*log(numurns)])
title('entropy as a function of spin')
xlabel('spin')
ylabel('S/Nk')