One Dimensional Diffusion Model

The model we would like to study today is the transport or diffusion of a material in a one-dimensional medium. This type of model is important in a variety of different physical situations: some examples include diffusion of a chemical compound or dye through a linear medium (such as in liquid chromatography), the flow of electrons through a wire, or the transport of pollution in a river. The way that we will setup the problem today is to study the flow of water in a river basin as spring melting occurs, which is obviously an important topic during April 1997 in Minnesota.

Consider a model where there is some initial concentration (such as a dye, pollution, or in our case runoff water) located at some location. We will break up the medium (in this case, the river) into a series of cells - in our case we might consider each cell to be a 5 mile long section of a river. For example, consider a system with 6 different cells with an initial concentration of 100 at location 3, which we represent as V(3)=100

          -----------------------------------------------------
         |  V(1)  |  V(2)  |  V(3)  |  V(4)  |  V(5)  |  V(6)  |
          -----------------------------------------------------

After a time interval (maybe 1 day in our river model) the material may either remain at the same location or diffuse to the left or right with some probabilities. In an undriven diffusion model (nothing driving the flow to the left or right), there would be equal amounts (100/3) in locations 2, 3 and 4 after one time step (in other words, V(2), V(3) and V(4) will all have the value 33.333). In the river model, there will be only a very small diffusion to the left (upstream), a relatively small probability of remaining at the same location, and a much larger probability of moving one cell to the right. After the next time interval, the process will repeat itself, with the material in each of these three cells diffusing to the left, right or remaining in the same location. We can make a Matlab script which will perform the diffusion and plot the graph each time. The algorithm will be as follows:

1. Initialize the number of cells to use and the number of time steps
2. Set all cells to zero except for the midpoint which contains the initial concentration
3. Perform a for loop (for each time step) which does the following
1. NewValues is replaced with the current Values times the fraction remaining in the same cell.
2. For each cell, the quantity in Values(j) will either be added to the cell to the left NewValues(j-1) or the right NewValues(j+1) after multiplying by the appropriate probabilities
3. After each time step, plot the concentration versus position

The program diff1.m below will perform the diffusion and plot the results after each iteration.

% DIFF1.M    1-Dimensional Diffusion Solution
%
% Tom Huber, March 1997
%
ncells = 25;       % Number of Cells in the Array
NTimes = 20;       % Number of Time steps to perform

FracLeft  = 1/3;   % Fraction diffusing to the left
FracSame  = 1/3;   % Fraction which remains in same cell
FracRight = 1/3;   % Fraction diffusing to the right

Values = zeros(ncells,1);  % Initialize the Array to have zero concentration
VMax = 100;                % Maximum Concentration
Values(ncells/2) = VMax;   % Initially set the concentration at the midpoint

for i=1:NTimes  % Perform a total of NTimes time steps
  NewValues = Values*FracSame;  % New values initially a fraction of original
  for j=2:ncells-1  % For each cell except leftmost or rightmost
        NewValues(j-1) = NewValues(j-1)+Values(j)*FracLeft;  % Diffuse Left
        NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
  end % for j=1...

  Values = NewValues;     % The updated values become the current values
  plot(Values)            % Plot the values
  axis([1 ncells 0 VMax]) % Set the minima and maxima on axes
  xlabel('Position')      % Label the axes and the title
  ylabel('Value')
  title(['After Time Step: ' num2str(i)])
  drawnow                % Make Matlab display the graph
                         % (Normally it only displays at the end of a program)
end  % for i=1...