Chapter 2
Iterative Methods for Solving Sets of Equations
2.1 The Gauss-Seidel Method
The Gauss-Seidel method may be used to solve a set of linear or nonlinear algebraic equations. We will illustrate the method by solving a heat transfer problem. For steady state, no heat generation, and constant k, the heat conduction equation is simplified to Laplace equation
2T = 0
For 2-dimensional heat transfer in Cartesian coordinate
+ = 0
The above equation can be put in the finite difference form. We divide the medium of interest into a number of small regions and apply the heat equation to these regions. Each sub-region is assigned a reference point called a node or a nodal point. The average temperature of a nodal point is then calculated by solving the resulting equations from the energy balance. Accurate solutions can be obtained by choosing a fine mesh with a large number of nodes. We will discuss an example from Incropera’s[1] text to illustrate the method.
Example 2.1-1
A long column with thermal conductivity k = 1 W/moK is maintained at 500oK on three surfaces while the remaining surface is exposed to a convective environment with h = 10 W/m2oK and fluid temperature T. The cross sectional area of the column is 1 m by 1 m. Using a grid spacing x = y = 0.25 m, determine the steady-state temperature distribution in the column and the heat flow to the fluid per unit length of the column.
Solution
The cross sectional area of the column is divided into many sub-areas called a grid or nodal network with 25 nodes as shown in Figure 2.1-1. There are 12 nodal points with unknown temperature, however only 8 unknowns need to be solved due to symmetry so that the nodes to the left of the centerline are the same as those to the right.
Figure 2.1-1 The grid for the column cross sectional area.
The energy balance is now applied to the control volume xy1 belongs to node 1 which is an interior node. To make the derivation general, node 1 can be considered as a node with index (m, n) in a two-dimensional grid as shown in Figure 2.1-1. The directions of conduction heat flow are assumed to be the positive x and y directions. For steady state with no heat generation
q(m-1, n)(m, n) + q(m, n-1)(m, n) = q(m, n)(m+1, n) + q(m, n)(m, n+1) (2.1-1)
where q(m-1, n)(m, n) is the conduction heat flow between nodes (m-1, n) and (m, n). Fourier’s law can be used to obtain
q(m-1, n)(m, n) = k(y1)
where (y1) is the heat transfer area with a unit depth and is the finite-difference approximation to the temperature gradient at the boundary between the two nodes. Appling Fourier’s law to each term in Equation (2.1-1) yields
k(y1) + k(x1) = k(y1) + k(x1)
The equation is divided by k(xy1) and simplified to
+ = 0(2.1-2)
For x = y, Eq. (2.1-2) becomes
Tm,n = (2.1-3)
The above result shows that the temperature of an interior node is just the average of the temperatures of the four adjoining nodal points. Using this formula, the temperatures for the first six nodes are
T1 = ( T2 + T3 + 500 + 500)
T2 = ( T1 + T4 + T1 + 500)
T3 = ( T1 + T4 + T5 + 500)
T4 = ( T2 + T3 + T6 + T3)
T5 = ( T3 + T6 + T7 + 500)
T6 = ( T4 + T5 + T8 + T5)
Nodes 7 and 8 are not interior points, therefore Eq. (2.1-3) is not applicable.
Figure 2.1-2 Directions of heat flow for nodes 7 and 8.
Making energy balance for node 7 yields
k(1) + k(x1) = k(1) + h(x1)(T7 300)
Multiplying the above equation by we obtain
500 T7 + 2T5 2T7 = T7 T8 + (T7 300)
= = 5.0
T7 = (2T5 + T8 + 2000)
Similarly an energy balance for node 8 yields
k(1) + k(x1) = k(1) + h(x1)(T8 300)
Multiplying the above equation by we obtain
T7 T8 + 2T6 2T8 = T8 T7 + (T8 300)
T8 = (2T6 + 2T7 + 1500)
We have 8 linear equations with 8 unknowns that can be solved by matrix method or iterations. Table 2.1-1 lists the Matlab program using Gauss-Seidel iteration to solve for the temperatures.
Table 2.1-1 ______
% Finite
% Initial guesses for the temperatures
T=400*ones(8,1);
for i=1:100
Tsave=T;
T(1)=.25*(T(2)+T(3)+1000);
T(2)=.25*(2*T(1)+T(4)+500);
T(3)=.25*(T(1)+T(4)+T(5)+500);
T(4)=.25*(T(2)+2*T(3)+T(6));
T(5)=.25*(T(3)+T(6)+T(7)+500);
T(6)=.25*(T(4)+2*T(5)+T(8));
T(7)=(2*T(5)+T(8)+2000)/9;
T(8)=(2*T(6)+2*T(7)+1500)/9;
eT=abs(T-Tsave);
if max(eT)<.01, break, end
end
fprintf('# of iteration = %g\n',i)
for i=1:8
fprintf('Node %g, Temperature = %8.2f\n',i,T(i))
end
> finite
# of iteration = 13
Node 1, Temperature = 489.30
Node 2, Temperature = 485.15
Node 3, Temperature = 472.06
Node 4, Temperature = 462.00
Node 5, Temperature = 436.95
Node 6, Temperature = 418.73
Node 7, Temperature = 356.99
Node 8, Temperature = 339.05
The heat transfer rate per unit depth from the column to the fluid is given as
q’/L = 2h[(500 300) + x(T7 300) + (T8 300)]
q’/L = (2)(10)(0.25)[100 + (356.99 300) + (0.5)(339.05 300)]
q’/L =883 W/m
If the matrix is strictly diagonally dominant, the Gauss-Seidel method will converge. However, convergence may also be achieved in many situations for which diagonal dominance cannot be obtained, although the rate of convergence is slowed. An nn matrix A is diagonally dominant if and only if
|ai,i| > ,i = 1, 2, , n
An example of a diagonally dominant system is the following set of equations
6x1 / 2x2 / + x3 / = / 11 / (6 > 2 +1)2x1 / + 7x2 / +2x3 / = / 5 / (7 > 2 +2)
x1 / + 2x2 / 5x3 / = / 1 / (5 > 1 + 2)
If possible, the equations should be reordered to provide diagonal elements whose magnitudes are larger than those of other elements in the same row. The Gauss-Seidel method may also be used to solve a set of nonlinear equations as shown in the next example.
Example 2.1-2
Use Gauss-Seidel method with the initial guess x = [0.1 0.1 0.1] to obtain the solutions to the following equations[2]
3x1 cos(x2 x3) = 0
81(x2 + 0.1)2 + sin x3 + 1.06 = 0
+ 20x3 + = 0
Solution
The above equations can be rearranged as follows
x1 = cos(x2 x3) +
(x2 + 0.1)2 = ( + sin x3 + 1.06) x2 = ( + sin x3 + 1.06)1/2 0.1
x3 =
Table 2.1-2 lists the Matlab program and the results using Gauss-Seidel iteration to obtain the solutions.
Table 2.1-2 ______
% Example 2.1-2 (Ex2d1d2), Gauss-Seidel iteration
%
x=[.1 .1 -.1];
for i=1:100
xsave=x;
x(1)=cos(x(2)*x(3))/3+1/6;
x(2)=sqrt(x(1)*x(1)+sin(x(3))+1.06)/9-0.1;
x(3)=-exp(-x(1)*x(2))/20-(10*pi-3)/60;
if abs(max(xsave-x))<1e-5, break, end
fprintf('x = %11.5f %11.5f %11.5f\n',x)
end
fprintf('# of iterations = %g\n',i)
fprintf('x = %11.5f %11.5f %11.5f\n',x)
> ex2d1d2
x = 0.49998 0.02223 -0.52305
x = 0.49998 0.00003 -0.52360
x = 0.50000 0.00000 -0.52360
# of iterations = 4
x = 0.50000 0.00000 -0.52360
1
[1] Incropera, F. P. and DeWitt, D. P., Fundamentals of Heat and Mass Transfer, John Wiley & Son, 2002.
[2]Numerical Analysis by Burden and Faires