## Question

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...

Tutorial Tasks

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.

## Solution Preview

This material may consist of step-by-step explanations on how to solve a problem or examples of proper writing, including the use of citations, references, bibliographies, and formatting. This material is made available for the sole purpose of studying and learning - misuse is strictly forbidden.

ncells = 100; % Number of Cells in the ArrayNTimes = 100; % Number of Time steps to perform

FracLeft = 1/10; % Fraction diffusing to the left

FracSame = 3/10; % Fraction which remains in same cell

FracRight = 6/10; % Fraction diffusing to the right

Values = zeros(ncells,1); % Initialize the Array to have zero concentration

VMax = 100; % Maximum Concentration

Values(ceil(ncells/4)) = 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......