############################################################

# Equilibrium temperature distribution
# (Laplace's equation and the mean value principle)

# Jim Carlson 
# http://www.math.utah.edu/~carlson)
# December 1, 1994

# THE PROBLEM:

# Find the interior temperatures in a thin metal plate
# at thermal equilibrium given the temperatures around 
# the edge of the plate.  One does this by breaking up 
# the plate into small square cells and using the mean value 
# principle: the temperature at a given cell is equal to the
# average temperature of its four nearest neighbors.

# We will try a 10 by 10 grid of cells where the temperature
# is 1 on two adjacent edges and 0 on the other two.


############################################################

  n := 10;  # number of cells in a row or column
  t := array(1..n, 1..n):  # array of temperatures

############################################################

# Set up an array of equations which express the mean value
# principle.

  eq := array(2..n-1, 2..n-1);  

  for i from 2 to n-1 do
    for j from 2 to n-1 do
      eq[i,j] := t[i,j] = (1/4)*( t[i-1,j] + t[i,j-1] + t[i+1,j] + t[i,j+1] );
    od;
  od;

############################################################

#  Define the temperature on the edge of the plate
#  (boundary conditions)

   # top and left
    for i from 1 to n
    do
      t[1,i] := 1.0;
      t[i,1] := 1.0;
    od:

    # bottom and right
    for i from 2 to n
    do
      t[n,i] := 0.0;
      t[i,n] := 0.0;
   od:

############################################################

#  Solve the equations and store the resulting 
#  temperatures in the array t.  We have to
#  convert the array of equations to a set
#  of equations because that is what solve needs.
#  We use the assign function to set the values
#  of the temperature array t.

   eqs := convert(eq, set):  
   sol := solve( eqs ):
   assign( sol );

############################################################

# Load the plots package and plot the temperatures.
# Suggestion: set the plot option to "patch and contour",
# set the coloration to "color by z value", and set
# the viewing position so that you look straight down
# the z axis.  This gives you a contour map where
# temperature is coded by color (red = hot, cold = blue,
# of course).

   with( plots );      # load package
   matrixplot( t );

############################################################