Temperature at Equilibrium --- The Discrete Heat Equation


Consider the problem of determining the temperature at interior points of a thin square plate given the temperature along the edges, assuming that the system is at equilibrium.  Starting from a quite simple and intuitive idea, this problem can be formulated in a way that leads immediately to a system of equations that can be solved using Gaussian elimination.  The key piece of physics is the mean value principle:
 

The temperature at equilibrium of a given point in a body at equilibrium is equal to the average temperature of nearby points.

 

 

To analyze our problem in an elementary way, we discretize it by dividing the square into a grid of many small "cells" as in the figure below:
 

          ________________________________
                 |       |       |       
            1    |   1   |   1   |   1
                 |       |       | 
          _______|_______|_______|________
                 |       |       |       
            1    |   a   |   b   |   0
                 |       |       |
          _______|_______|_______|________
                 |       |       |       
            1    |   c   |   d   |   0
                 |       |       |
          _______|_______|_______|________
                 |       |       |       
            1    |   0   |   0   |   0
                 |       |       |
          _______|_______|_______|________
We have labeled each cell by a number or letter.  The numbers around the edge of the square plate represent measured temperatures.  The letters in the interior cells represent unmeasured temperatures which we are to solve for.  Notice that it is hot on the top and left of the square plate and cold on the bottom and the right.

We can deduce the interior temperatures from the measured boundary temperatures since the mean value principle tells us that
 

The temperature of an interior cell is equal to one fourth of the sum of its four nearest neighbors.

This gives four equations in the four interior temperatures a,b,c and d:
               a = (1 + 1 + b + c)/4
               b = (a + 1 + 0 + d)/4
               c = (1 + a + d + 0)/4
               d = (c + b + 0 + 0)/4
Thus we can find the temperatures by row reducing the augmented matrix corresponding to this system.
 
     >> M = [   1   -1/4   -1/4     0    1/2;
              -1/4    1      0    -1/4   1/4;
              -1/4    0      1    -1/4   1/4;
                0   -1/4   -1/4     1     0   ]
     >> R = rref(M)
At this point the last column of the reduced matrix R holds the interior temperatures we computed.   We find that a = 3/4, b = 1/2, c = 1/2 and d = 1/4 at equilibrium.     Next we illustrate one method for visualizing the temperature distribution since numerical data is often easier to understand if we can visualize it.  We want to set up a 4 x 4 matrix T that holds the 16 temperatures for our 4 x 4 grid.   One way to do this is to set up a zero matrix of the correct shape and then modify its entries as follows:
 
     >> T = zeros(4,4)
     >> T(:,1) = 1               % set all entries in column 1 to be 1
     >> T(1,:) = 1               % set all entries in row 1 to be 1
     >> T(2:3,2:3) = reshape( R(:,5), 2, 2 )'


This last command needs explanation.  The MATLAB function on the right-hand-side reads down column 5 of the reduced matrix R and reshapes it as a 2 x 2 matrix.  It makes the first two entries into the first column of the new 2 x 2 matrix and the next two entries become the second column.    We actually want the first two entries of column 5 to be the first row of the 2 x 2 matrix T, so we should take the transpose.  In MATLAB we use the single quote to transpose a matrix.  (In this example it does not matter if we transpose - with our choice of boundary condition the solution will automatically be symmetric.)  For more information about the reshape command, go to the helpdesk and look up this MATLAB function there.

We can get  a plot of the temperatures using the MATLAB command

     >> imagesc(T);
Here color represents temperature, with hotter regions represented with more red and colder regions represented with more blue tones.

Of course, our results with a 4 x 4 grid are very crude.   A finer grid will give us a better approximation to the equilibrium temperature distribution, but of course we have to work a bit harder to do the computation.  For example, if we use a 5 x 5 grid instead of a 4 x 4 we go from 4 equations in 4 variables to 9 equations in 9 variables:
 

          ________________________________
               |      |      |      |
            1  |   1  |   1  |   1  |  1
          _____|______|______|______|_____
               |      |      |      | 
            1  |   a  |   b  |   c  |  0
          _____|______|______|______|_____
               |      |      |      | 
            1  |   d  |   e  |   f  |  0
          _____|______|______|______|_____
               |      |      |      |  
            1  |   g  |   h  |   i  |  0
          _____|______|______|______|_____
               |      |      |      |  
            1  |   0  |   0  |   0  |  0
          _____|______|______|______|_____


Exercise 1:  Set up the system of 9 equations in 9 variables corresponding to the 5 x 5 grid pictured above.  Solve for the 9 interior temperatures.  Use the colon operator and the reshape command to define a matrix T that represents the 25 temperatures as illustrated above and plot with imagesc.  (Of course, you must modify the sample MATLAB commands above to reflect the fact that the augmented matrix M will now be a 9 x 10 matrix, R should be its reduced row echelon form,  and now T will be a 5 x 5 matrix of temperatures.)
 

Larger Grids and a Little Programming

When we discretize our temperature problem by dividing  the square into cells we get an approximate temperature distribution, not the exact one.  The more cells we use, the finer the subdivision, the more accurate the approximation.  Suppose for better accuracy we use a mesh half as big as the 5 x 5 grid.  This means a 10 x 10 grid and thus 64 equations in 64 unknowns.  Without a computer it is not feasible to write out such big systems, much less solve them.  And even with a solution in hand, a list of the 64 interior temperatures, we would have a much easier time understanding it if we could easily view the data in graphical form.

Nonetheless, it is possible to work with very large systems.  The key point is that the mean value equations follow a regular pattern and so can be defined by a short MATLAB program.  To see how this is done, recall once again that the temperature of an interior cell at equilibrium is the average temperature of its four nearest neighbors.  These are the cell directly above, the cell directly below,  the next cell to the left and the next cell to right.  We will use an n x n matrix T to represent the temperatures.  If we consider the temperature T(i,j) in the cell in row i and column j of our grid, the cell directly above lies in row i-1 and column j.  The cell directly to the right lies in row i and column j + 1 and so on.  Thus the interior temperatures must satisfy equations of the form

        T(i,j) = ( T(i-1,j) + T(i,j+1) + T(i+1,j) + T(i,j-1) )/ 4
We start by setting up a matrix T that will eventually hold all the temperature values, but to get started we will define a matrix that has the correct temperatures along the boundary and zeros everywhere else.
     >> n = 5;                  % we'll do the 5 x 5 case as an example
     >> vars = (n-2)^2;         % an n x n grid gives (n-2)^2 variables
     >> T = zeros(n)            % set up an n x n zero matrix
     >> T(1,:) = 1.0            % set all temperatures in row 1 to be 1.0
     >> T(:,1) = 1.0            % set all temperatures in column 1 to be 1.0
At this point, T holds the correct boundary temperatures but we still need to compute the nine unknown interior temperatures
                 T(2,2)     T(2,3)     T(2,4)
                 T(3,2)     T(3,3)     T(3,4)
                 T(4,2)     T(4,3)     T(4,4)
We can loop through these interior cells using a double loop in MATLAB.
     >> for row = 2:n-1
           for col = 2:n-1
               -------  need to fill in some code here  -----------
               -------  to set up eqn corresponding to  -----------
               -------      each interior cell          -----------
               -------    see explanations below        -----------
           end
        end
It will be easier to think of these 9 unknown temperatures as a single vector.  The ordering on this vector should match the order we go through the cells in the double loop above.
                 T(2,2) = x_1      T(2,3) = x_2      T(2,4) = x_3
                 T(3,2) = x_4      T(3,3) = x_5      T(3,4) = x_6
                 T(4,2) = x_7      T(4,3) = x_8      T(4,4) = x_9
Our problem is to set up the augmented matrix  M that represents the corresponding system of nine equations in these nine unknown temperatures.  The equation corresponding to the 2,2 cell tells us
           T(2,2) = ( T(1,2) + T(2,3) + T(3,2) + T(2,1) )/4
or in other words, taking the boundary temperatures T(1,2) = 1 and T(2,1) = 1 into account, we have
               x_1  - (1/4)x_2  - (1/4)x_4 = 1/2
In general, for the cell corresponding to x_k we will get an equation where x_k appears with coefficient 1.  The next cell to the right will either be a boundary cell or it corresponds to the variable x_(k+1).  If it is a boundary cell then we adjust the constant term and otherwise we know that x_(k+1) appears in this equation with coefficient -1/4.  The cells directly above and below are a little more complicated.  In the 5 x 5 grid case, the cell directly above will either be a boundary cell or it corresponds to the variable x_(k-3).  More generally,  the n x n grid gives rise to (n-2)^2 equations in (n-2)^2 variables and the cell directly above the kth cell will in this case correspond to x_(k - (n-2)).  With these ideas we can automate the process to set up the augmented matrix M:

 We start by setting up a zero matrix of the correct shape, in the 5 x 5 case it will be a 9 x 10 matrix:

     >> M = zeros( vars, vars + 1)
Next we will loop through all the interior cells, each of which corresponds to a variable x_k.  The mean value principle gives us a linear equation involving x_k and perhaps the variables that correspond to the four adjacent cells.   For each value of k we must consider the four adjacent cells separately.     If  it is a boundary cell we adjust the constant term (that is, the last column of M) and otherwise, we set the coefficient corresponding to that interior cell to be -1/4:
     >>  k = 0;
         for row = 2:n-1
            for col = 2:n-1
                 k = k + 1;             % keep track of which equation we're working on
                 M(k,k) = 1             % since x_k occurs with coefficient 1 in the kth eqn
                 % Now deal with the cell directly above the kth cell.  How does it contribute?
                 % It is a boundary cell if we are in the top interior row, that is if row = 2.
                 % In that case, adjust the constant term, stored in the last col of M.
                 % Otherwise, the cell directly above corresponds to the variable x_(k-n+2)
                 % and it appears with coefficient - 1/4
                 if row == 2
                   M(k, vars + 1) = M(k,vars + 1) + T(row-1,col)/4;
                 else
                    M(k, k-n+2) = -1/4;
                 end;
                 % Now you should be able to add the code to deal with the cell directly
                 % below, the next cell to the left and the next cell to the right by making
                 % small modifications to the if-else-end statement above.
           end;
         end;
Exercise 2:   Complete the MATLAB code above to automate the solution using a 5x5 grid.  Use your work in Exercise 1 as a check.

Exercise 3:  Rework the problem with a 10 x 10 grid.  What do you expect to happen as we refine the grid?  Can you predict the limit?  (You might want to take a look at the graph produced by imagesc, but don't print it or turn it in.)