% VALSURF.M Plot the surface from the diffusion model
clf; % Clear the current figure
surf(AllValues) % Plot the surface
colormap(jet) % Change the colormap
view(0,90) % Set the viewpoint to be above
colorbar % Show a colorbar which has the values
xlabel('Time Step') % Label the axes
ylabel('Position')
% plotloc.m Plot Results at some location
% Uses the array AllValues and plots the results
%
loc = input('Enter Cell Location to Display ');
Values = AllValues(loc,:);
plot(Values)
xlabel('Time Step')
ylabel('Concentration')
title(['Cell Number ' num2str(loc)])
[ValMax LocMax] = max(Values);
disp(['Max Value = ' num2str(ValMax) ' At Time Step ' num2str(LocMax)])
% plotloc2.m Plot Results at some location
% Uses the array AllValues and plots the results
%
loc = input('Enter Cell Location to Display ');
clf
Value1 = AllValue1(loc,:);
plot(Value1,'g-')
Values = AllValues(loc,:);
hold on
plot(Values,'r-')
xlabel('Time Step')
ylabel('Concentration')
title(['Cell Number ' num2str(loc)])
legend('Model 1','Model 2')
[ValMax LocMax] = max(Value1);
disp(['Model 1 Max Value = ' num2str(ValMax) ...
' At Time Step ' num2str(LocMax)])
[ValMax LocMax] = max(Values);
disp(['Model 2 Max Value = ' num2str(ValMax) ...
' At Time Step ' num2str(LocMax)])
% 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 Y-Data 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 1-Dimensional 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 = (1-MeltFrac)*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(j-1) = NewValues(j-1)+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 1-Dimensional 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 4-8) % with 20 instead of 1 reservoir with 100 VRemain(4:8) = 20*ones(5,1); Values = MeltFrac*VRemain; % Initially let some fraction melt VRemain = (1-MeltFrac)*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(j-1) = NewValues(j-1)+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
% DIFF1VEC.M 1-Dimensional 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 4-8) % with 20 instead of 1 reservoir with 100 VRemain(4:8) = 20*ones(5,1); Values = MeltFrac*VRemain; % Initially let some fraction melt VRemain = (1-MeltFrac)*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 ncells-1 and % add them either to the left or right NewValues(1:ncells-1)=NewValues(1:ncells-1)+FracLeft(2:ncells).*Values(2:ncells); NewValues(2:ncells)=NewValues(2:ncells)+FracRight(1:ncells-1).*Values(1:ncells-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...
Electronic Copy: http://physics.gac.edu/~huber/envision/tutor2/diffvian.htm
Created: 4-APR-1997 by Tom Huber, Physics Department, Gustavus Adolphus College.