To view the results as a movie, you can use the
file Valmovie.m below % Valmovie.m Plot Results for successive times
[NCells NSteps] = size(AllValues);
VMax = max(max(AllValues));
clf
for i=1:NSteps
Values = AllValues(:,i);
plot(Values,'r')
xlabel('Location')
ylabel('Concentration')
title(['Time Step ' num2str(i)])
axis([1 NCells 0 VMax])
drawnow
end
Another variation on Valmovie.m (which I call fstmovie.m
uses Matlab's ability to do animation using HandleGraphics.
In this program, instead of displaying the axes & titles each time,
only the YData values are changed on the graph. The resulting
movie executes much faster than the previous version
% Fstmovie.m Plot Results for successive times
% Uses HandleGraphics for faster animation of the movie
[NCells NSteps] = size(AllValues);
VMax = max(max(AllValues));
TDelay = .1; % This is the minimum amount of time between movies
clf;
plot ([],[]);
axis([1 npts 0 VMax]);
axis(AXIS);
x = 1:NCells; % Initialize X Data for 1 to NCells
% Draw the line and store its handle in LineData
LineData = line('color','y','erase','xor','xdata',x,'ydata',x);
xlabel('Position')
ylabel('Concentration')
% Display Iteration number and store its handle in TextData
TextData = text(npts/10,0.9*VMax,'Iteration: 0','Color','r');
set(TextData,'erase','Background')
drawnow
for i=1:NSteps
Time0 = clock; % Get the current time
Values = AllValues(:,i);
set(LineData,'ydata',Values); % Update the Y Values on the Graph
set(TextData,'string',['Iteration: ' num2str(i)]) % Update text
drawnow
while etime(clock,Time0)<TDelay % Pauses at least .1 second
% disp('This pauses on a fast CPU')
end % while ...
end % for i...
Group Tutorial Tasks

The first task has a reservoir that "melts" at some rate (this is
similar to the model melting.m which we did
at the beginning of the Tutorial on Matlab Scripts.

Comparing the runs shown here, one notes that a slower melting
rate leads to a lower peak value, but also leads to a longer flooding
time (the peak is wider)

% DIFF1MLT.M 1Dimensional Diffusion Solution with a Melting Reservoir
%
% Tom Huber, March 1997
%
ncells = 25; % Number of Cells in the Array
NTimes = 35; % 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
MeltFrac = 0.1; % Fraction which "Melts" each day
VRemain = VMax;
Values(ncells/4) = MeltFrac*VRemain; % Initially let some fraction melt
VRemain = (1MeltFrac)*VRemain; % This is fraction remaining after step
AllValues = Values;
for i=1:NTimes % Perform a total of NTimes time steps
% First we let the current water flow downstrem
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(j1) = NewValues(j1)+Values(j)*FracLeft; % Diffuse Left
end
if j<ncells
NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
end
end % for j=1...
Values = NewValues; % The updated values become the current values
% Now we will add some new melted water to the river
Melted = MeltFrac*VRemain; % This is how much melted
VRemain = VRemain  Melted; % This is how much remains
Values(ncells/4) = Values(ncells/4) + Melted; % Add the melted
AllValues = [AllValues Values];
end % for i=1...
 The next task is to have a reservoir at every location.
This is accomplished
by making VRemain and Meleted vectors.
In the example below, the reservoir of 100 is replaced by 5 reservoirs of
20 (the melting fraction is 0.4 in both cases).

Comparing the runs shown here, one notes that an extended
reservoir (5 cells of 20 each  solid curve) leads to a
lower peak value, but also leads to a longer flooding
time than a single reservoir

% DIFF1RES.M 1Dimensional Diffusion Solution with a Reservoir
% at each location
%
% Tom Huber, March 1997
%
ncells = 25; % Number of Cells in the Array
NTimes = 35; % 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
MeltFrac = 0.4; % Fraction which "Melts" each day
VRemain = zeros(ncells,1); % Initialize all reservoirs to zero
% Now load 5 reservoirs (locations 48)
% with 20 instead of 1 reservoir with 100
VRemain(4:8) = 20*ones(5,1);
Values = MeltFrac*VRemain; % Initially let some fraction melt
VRemain = (1MeltFrac)*VRemain; % This is fraction remaining after step
AllValues = Values;
for i=1:NTimes % Perform a total of NTimes time steps
% First we let the current water flow downstrem
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(j1) = NewValues(j1)+Values(j)*FracLeft; % Diffuse Left
end
if j<ncells
NewValues(j+1) = NewValues(j+1)+Values(j)*FracRight; % Diffuse Right
end
end % for j=1...
Values = NewValues; % The updated values become the current values
% Now we will add some new melted water to the river
Melted = MeltFrac*VRemain; % This is how much melted
VRemain = VRemain  Melted; % This is how much remains
Values = Values + Melted; % Add the melted
AllValues = [AllValues Values];
end % for i=1...
 Finally, the left, right and same fractions can be treated
as arrays and a constriction can be placed some location (in this case cell
10). After running diff1vec.m,
view the results using valsurf and valmovie

Comparing the runs shown here, one notes that the crest
value and the length of flooding is significantly higher at the
location of the constriction (in this case the constriction and the
plot are both at location 10).


Comparing the runs shown here, one notes that the crest
value is significantly lower and occurs later
downstream of the constriction (in this case the constriction is
at location 10 and the
plots is at location 10).


The surface plot (from valsurf) also allows one to see
the effect of the constriction be noticing that diffusion
occurs very slowly beyond cell 10

% DIFF1VEC.M 1Dimensional Diffusion Solution with a Reservoir
% at each location and the diffusion fractions
% stored in vectors
%
% Tom Huber, March 1997
%
ncells = 25; % Number of Cells in the Array
NTimes = 35; % Number of Time steps to perform
FracLeft = 1/10*ones(ncells,1); % Fraction diffusing to the left
FracSame = 3/10*ones(ncells,1); % Fraction which remains in same cell
FracRight = 6/10*ones(ncells,1); % Fraction diffusing to the right
% Make a constriction at one cell (cell number 10)
FracLeft(10) = 2/10;
FracSame(10) = 6/10;
FracRight(10) = 2/10;
MeltFrac = 0.4; % Fraction which "Melts" each day
VRemain = zeros(ncells,1); % Initialize all reservoirs to zero
% Now load 5 reservoirs (locations 48)
% with 20 instead of 1 reservoir with 100
VRemain(4:8) = 20*ones(5,1);
Values = MeltFrac*VRemain; % Initially let some fraction melt
VRemain = (1MeltFrac)*VRemain; % This is fraction remaining after step
AllValues = Values;
for i=1:NTimes % Perform a total of NTimes time steps
% First we let the current water flow downstrem
% NOTE: Sinc FracSame is now a vector, we use the .* operator
NewValues = Values.*FracSame; % New values initially a fraction of original
% NOTE: We must use subsets of the vectors of length ncells1 and
% add them either to the left or right
NewValues(1:ncells1)=NewValues(1:ncells1)+FracLeft(2:ncells).*Values(2:ncells);
NewValues(2:ncells)=NewValues(2:ncells)+FracRight(1:ncells1).*Values(1:ncells1);
Values = NewValues; % The updated values become the current values
% Now we will add some new melted water to the river
Melted = MeltFrac*VRemain; % This is how much melted
VRemain = VRemain  Melted; % This is how much remains
Values = Values + Melted; % Add the melted
AllValues = [AllValues Values];
end % for i=1...