Answers - Visualization of Diffusion Model Results

In the April 1997 tutorial there were many tasks for you to perform. This document shows some script files and plots for solving these tasks. It is important to note that there are often many different ways to solve the problems - it is likely that your script files will use a different approach

To run these script files, you must save them as Matlab m-files. One way is to highlight the section of code, copy it to the clipboard, paste it into a new document and save it as a file (such as diff1b.m). The other approach would be to save the entire file diffvian.htm as a file and then edit out each script file. The only problem with this approach is in dealing with the less-than (<) and greater-than (>) symbols. Because of how HTML interprets these symbols, I have replaced all < and > symbols with the text strings &lt; and &gt;. You will have to edit the source code and replace the string &lt; with < (and &gt; with >)

Copy your program program diff1.m to diff1b.m and modify it as shown below.
      % DIFF1B.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
      AllValues = Values;

      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<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
        AllValues = [AllValues Values];
      end  % for i=1...

Tutorial Tasks

  1. Make a script file valsurf.m that will plot and label the surface
              % 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')
    
  2. Extracting regions of the data set Values=AllValues(10,:) and Values=AllValues(:,5) the results are shown below
    AllValues(10,:) Value = AllValues(10,:)
    This can be interpreted as the concentration at the 10th location. (Click on image for an enlarged view)
    AllValues(:,5) Value = AllValues(:,5)
    This can be interpreted as the concentration after the 5th time step
  3. Make an M-File plotloc.m to plot the results at a selected location
          % 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)])
  4. When does the crest (maximum) occur for a position 5 cells downstream from the initial concentration? What about 15 cells downstream? For 5 cells downstream from maximum (at Location 10), the crest value is 22.6 and occurs at time step 7. For 15 cells downstream, we need to extend the time range in diff1b.m by setting NTimes=35. This results in a crest value of 11.4 at time step 27. Does doubling the initial concentration double the crest at these locations? Yes
  5. Running the model with different left, right, and same fractions, we can view the results using plotloc2.m below.
    AllValues 2 Plots Comparing the original run (Model 1) with left, same, right fractions of (1/10, 3/10, 6/10) to a run (Model 2) with fractions (0.5/10, 3/10, 6.5/10) we see that the peak in the second model comes earlier and is higher.
          % 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)])
    
    
  6. 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 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

    1. 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.
      2 Melting Plots 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...
    2. 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).
      Extended Reservoir 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...
      
    3. 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
      Location of Constriction 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).
      Downstream of Constriction 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).
      Surface Plot 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    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.