interface( warnlevel=0 ):
with( laylinalg ):
# Text Reference: Section 2.5 and Exercise 2.5-31
# Lab 4: Equilibrium Temperature Distributions
# The purpose of this Lab is to discuss the equilibrium temperature of a thin plate, a problem which leads to a system of linear equations. Methods for solving these equations will be compared, and some methods will be shown to be more efficient than others.
#
# Consider a thin square plate whose faces are insulated from heat. Suppose that the temperature along the four edges of the plate is known, and further suppose that those temperatures are held constant. After some time has passed, the temperature inside the plate will reach an equilibrium. Finding this equilibrium temperature distribution at the points on the plate is desirable, given only the temperature data from the edges of the plate. Unfortunately, the exact determination of this temperature distribution is a difficult problem. An approximation to the exact distribution may be found by discretizing the problem; that is, by only considering a few points on the plate and approximating the temperature at those points.
#
# A property from thermodynamics helps with the discretization of the problem. The property says that if equilibrium has been achieved, then the temperature at a point is the average value of the temperature at surrounding points.
# Sample modeling for a rectangular dam and corresponding calculations with the Mean Value Property can be found here:
# http://www.math.utah.edu/~gustafso/s2015/2250/quiz/sampleQuizzes/quiz7/quiz7-sample.pdf
# The Mean Value Property:
# If a plate has reached thermal equilibrium, and P is a point on the plate, and C is a circle centered at P fully contained in the plate, then the temperature at P is the average value of the temperature function over C. See Figure 1 for a picture of this situation.
#
# Figure 1
#
# In order to discretize the problem, place a grid over the plate and concentrate only on the points where the grid lines cross; the temperature at only those points will be considered. The grid is fashioned so that some grid points lie on the boundary of the plate; assume that the temperature at these points equals the external temperature. At grid points inside the plate, assume the following numerical version of the Mean Value Property.
#
# The Discretized Mean Value Property:
# If a plate has reached thermal equilibrium, and P is a grid point not on the boundary of the plate, then the temperature at P is the average of the temperatures at the four closest grid points to P. Assumed is an equally-spaced grid of points not on the boundary of the plate.
# Example:
# Consider placing the following grid on the square plate; see Figure 2. There are four points inside the plate to consider; the temperatures at these points are labelled x[1];, x[2];, x[3];, and x[4];. Assume that the exterior temperatures are as labelled in Figure 2.
#
#
# Figure 2
#
# The Discretized Mean Value Property gives rise to the following four equations. The temperatures of the four closest grid points are recorded clockwise in map order: West, North, East, South. For example, x[1];equals the the sum of four temperatures divided by 4 (the average of 4 temperatures), where West=10, North=30, East=x[2], ;South=x[3].;
# x[1] = 1/(4)(10+30+x[2] +x[3]),;
# x[2] = 1/(4)(x[1]+30+0+x[4] ),;
# x[3] = 1/(4)(10+ x[1] +x[4]+20),;
# x[4] = 1/(4)( x[3] +x[2]+0+20).;
# These four equations are equivalent to the following system of linear equations:
# 4 x[1] = x[2] +x[3]+40,;
# 4 x[2] = x[1] +x[4]+30,;
# 4 x[3] = x[1] +x[4]+30,;
# 4*x[4] = x[2]+x[3]+20.;
# What if the external temperatures change? A new system of equations would need to be solved, but it would be very similar to the one previously considered. To handle this difficulty more easily, the system may be rewritten in matrix-vector form, as we now show.
# Define x = Vector(4, {(1) = x[1], (2) = x[2], (3) = x[3], (4) = x[4]});, C = Matrix(4, 4, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 1, (1, 4) = 0, (2, 1) = 1, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (3, 1) = 1, (3, 2) = 0, (3, 3) = 0, (3, 4) = 1, (4, 1) = 0, (4, 2) = 1, (4, 3) = 1, (4, 4) = 0}); and b = Vector(4, {(1) = 40, (2) = 30, (3) = 30, (4) = 20}); , then 4x = Cx+b, and x is the vector of equilibrium temperatures. Notice that C functions as a type of adjacency matrix; the (i;,j;) entry is a 1 if the points corresponding to x[i]; and x[j]; are connected in the grid, and is 0 otherwise. Define A=4I-C and write 4x = Cx+b in the form Ax = b. The solution x to system Ax = b may be written as x=1/A;b.
#
# Even though this formula appears to simplify matters when the external temperatures change, difficulties with using 1/A; abound, especially for large systems. See the Numerical Notes in sections 2.2 and 2.5 (textbook 5/E, pages 111 and 129) for more details. The LU decomposition may be used to alleviate these problems, as in textbook exercise 2.5-31. If A=LU, the system Ax=b becomes LUx=b. , the systems Ly=x and Ux=y can be solved in turn to find x.
# Problems to be Submitted:
# Problem 1.
# Solve the above system Ax = b by row reduction (swap, scale, replace operations).
# Problem 2.
# Find A;, 1/A; and x for the above system Ax = b. Confirm that your answer matches that of Problem 1.
# Problem 3.
# Suppose that the external temperature at the top of the plate changes to 50°. Use results from Problem 2 to find the new equilibrium temperature distribution.
# Problem 4.
# Do an LU decomposition of A above. Use it to solve the original system and the system with altered external temperatures. Confirm your earlier results.
# Problem 5.
# A finer grid should give a better approximation of the equilibrium temperatures. Consider the grid in Figure 3.
#
# Figure 3
#
# Now there are 25 grid points inside the plate, so we now have a system of 25 equations in 25 unknowns. The Maple commands given below define the matrix C and the vector b for this problem. Find the temperature distribution in this case with the original external temperatures by either finding (4I - C)^(-1); or by using the LU factorization. Compare the refined answer to the answer from Problem 1 (e.g., grid point x[19];appears in Figure 1 as grid point x[4];).
#
# C := Matrix(
# [ [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1],
# [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0] ] );
#
# b := Vector( [40, 30, 30, 30, 30, 10, 0, 0, 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 0, 0, 30, 20, 20, 20, 20] );
#