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:

- Initialize the number of cells to use and the number of time steps
- Set all cells to zero except for the midpoint which contains the initial concentration
- Perform a
`for`loop (for each time step) which does the following `NewValues`is replaced with the current`Values`times the fraction remaining in the same cell.- 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 - 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...

- Run the M-file
`diff1`as written (this performs undriven diffusion, with equal probabilities to the left, right and same cells) - Modify the M-file to better model the flow in the river by starting the initial concentration a quarter of the way along the array and changing the fraction to the left to 1/10, the fraction staying in the same location to 3/10, and the fraction to the right to 6/10. Run the program several times with different values for the relative fractions, NMax, NCells, etc.
- Start the initial concentration in Cell number 1 - what
happens? The model as written has a problem with the
leftmost and rightmost points, namely we cannot add onto the points to
the left or to the right of the
endpoints [which would be
`NewValues(0)`or`NewValues(ncells+1)`respectively]. We got around this problem by making our loop`for j=2:ncells-1`, however this causes the program to miss diffusion to the right from the leftmost cell (and to the left from the rightmost cell). Correct this problem by making our loop run from`for j=1:ncells`and include`if`statements whereby we have transport to the left if`j>1`and transport to the right if`j<ncells.`

Created: 4-APR-1997 by