############################################################ # 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 ); ############################################################