% 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.`% 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/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(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=1:ncells % For each cell except leftmost or rightmost if j>1 NewValues(j-1) = NewValues(j-1)+Values(j)*FracLeft; % Diffuse Left end if j

% 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/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(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...

Created: 4-APR-1997 by Tom Huber, Physics Department, Gustavus Adolphus College.